Finite Field Polynomial Multiplication in the Frequency Domain with ...
Finite Field Polynomial Multiplication in the Frequency Domain with Application to Elliptic Curve Cryptography S. Baktır, B. Sunar June 12, 2005 Abstract We introduce an efficient method for computing Montgomery products of polynomials in the frequency domain. The discrete Fourier transform (DFT) based method originally proposed for integer multiplication provides an extremely efficient method with the best asymptotic complexity, i.e. O(m log m log log m), for multiplication of m-bit integers or (m−1)st degree polynomials. Unfortunately, the DFT method bears significant overhead due to the conversions between the time and the frequency domains. This makes the original DFT method impractical for multiplication of short operands as used in many applications. In this work, we introduce an algorithm which performs the entire modular multiplication (including the reduction step) in the frequency domain, and thus eliminates costly back and forth conversions between the frequency and the time domains. Furthermore, with careful selection of parameters, we show that in computational platforms where multiplication operation is expensive frequency domain multiplication of finite field elements may be realized more efficiently than multiplication in the time domain for operand sizes relevant to elliptic curve cryptography. Key Words: Finite fields, modular multiplication, discrete Fourier transform (DFT), number theoretic transform (NTT), elliptic curve cryptography (ECC).
1
Introduction
Finite fields have many applications in coding theory [5, 4], cryptography [15, 7, 17], and digital signal processing [6]. Hence efficient implementation of finite field arithmetic operations is crucial. Multiplication of field elements is commonly implemented in terms of modular multiplication of polynomials realized in two steps: ordinary polynomial multiplication and modular reduction of the result by the field generating polynomial. Unlike polynomial multiplication, modular reduction has only linear complexity, i.e. O(m), for polynomials of degree m − 1 when a fixed special modulus is chosen. Hence it is essential to attack the complexity of the multiplication step. The classical schoolbook polynomial multiplication method for multiplication of (m − 1)st degree polynomials has quadratic complexity, i.e. O(m2 ), given in terms of coefficient multiplications and additions. The complexity may be improved to O(mlog2 3 ) using the Karatsuba algorithm [11]. Despite the significant improvement gained by the Karatsuba algorithm, the complexity is still not optimal. Furthermore, the implementation of the Karatsuba algorithm is more burdensome due to its recursive nature. The known fastest multiplication algorithm, introduced by Sch¨ onhage and Strassen [21], performs multiplication in the frequency domain using the Fast Fourier Transform (FFT) [9] with complexity O(m log m log log m) for multiplication of m-bit integers or m-coefficient polynomials [10]. Unfortunately, the FFT based algorithm becomes efficient and useful in practice only for very large operands due to the overhead associated with the forward and inverse Fourier transform 1
operations. When the transformation computations are not considered, frequency domain polynomial multiplication (without modular reduction) has surprisingly low linear, i.e. O(m), complexity. Sadly, no efficient method for performing modular reduction is known to exist in the frequency domain, and therefore one needs to convert the product of a polynomial multiplication operation back to the time domain to perform modular reduction bearing significant overhead. Similar approaches on doing arithmetic over domains where multiplication has only linear complexity but reduction is difficult have been proposed with applications in cryptography [2],[3]. In [2], the authors propose a full implementation of RSA[20] in the Residue Number System (RNS) which does not require any conversions either from radix to RNS beforehand or from RNS to radix afterwards. In their implementation they use an RNS version of the Montgomery multiplication algorithm. In [3] an alternate polynomial representation is utilized where a polynomial is represented by its values at sufficiently many points and a Montgomery like modular multiplication algorithm is proposed for elliptic curve cryptography over the extension fields GF (pm ) for p > 2m. In this paper we propose a new algorithm for achieving modular reduction in the frequency domain. With our approach, the entire finite field multiplication, including modular reduction, can be carried out in the frequency domain. In many finite field applications a chain of arithmetic operations are performed rather than a solitary one. Using our method after an initial conversion step all intermediary operations may be computed in the frequency domain. Therefore, there will be no need for conversion before and after every single finite field multiplication except before the very first and after the very last ones. In the next section, we briefly give some background information on finite fields, the discrete Fourier transform (DFT) in finite fields and application of the DFT to multiplication. In Section 3, we propose a new algorithm, named DFT modular multiplication, for finite field modular multiplication of polynomials in the frequency domain. In Section 4, we present ideas for efficient implementation of DFT modular multiplication. Furthermore, in Section 5, we give a complexity analysis of the DFT modular multiplication algorithm, compare it with other efficient methods such as the Karatsuba algorithm, and show its relevance for elliptic curve cryptography. Finally, in Section 6 we present our conclusions.
2 2.1
Background Finite Fields and Polynomial Representation
A field with finite number of elements is called a finite field or Galois field, denoted by Fq or GF (q), where q stands for the number of elements in the field [14]. The number of elements in a finite field is always a prime or a prime power, i.e., q = p or q = pm ,where the prime number p is called the characteristic of the finite field. When q is a prime, i.e. q = p, the finite field GF (p) is called a prime field. The prime field GF (p) is the field of residue classes modulo p and its elements are represented by the integers in {0, 1, 2, . . . , p − 1}. When q is a prime power, i.e. q = pm , the finite field GF (pm ) is called an extension field. The extension field GF (pm ) is generated by using an mth degree irreducible polynomial over GF (p) and it is the field of residue classes modulo the irreducible field generating polynomial. Hence, in polynomial representation the elements of GF (pm ) are represented by polynomials of degree m − 1 with coefficients in GF (p).
2
2.2
Multiplication in Finite Fields Using the Schoolbook Method
Finite field polynomial multiplication is performed modulo Pm−1 thei field characteristic Pm−1and ithe field generating polynomial. For instance, let A(x) = a x and B(x) = i i=0 i=0 bi x , where ai , bi ∈ GF (q), be two polynomials representing two elements of the finite field GF (q m ) constructed with an mth degree field generating polynomial f (x) which is irreducible over GF (q). The most straightforward way of performing multiplication of two polynomials would be by the the schoolbook method of multiplication. One can perform the multiplication of A(x) and B(x) by the schoolbook method and find the product C(x) = A(x) · B(x) in two steps as follows: 1. Polynomial multiplication: 0
C (x) = A(x) · B(x) =
m−1 X m−1 X
ai bj xi+j
i=0 j=0
2. Modular reduction: C(x) = C 0 (x) (mod f (x))
Note above that C(x) is an (m − 1)st degree polynomial whose coefficients are elements of GF (q). As much as it is straightforward, schoolbook method of multiplication is inefficient and usually not the preferred method in finite field applications. In section 3, we will introduce an efficient method for finite field polynomial multiplication by using the Number Theoretic Transform [18], and compare it with other efficient algorithms such as the Karatsuba Multiplication [11] and a more recent method by Bajard et al. that utilizes an alternate polynomial representation and Lagrange’s theorem [3].
2.3
Number Theoretic Transform (NTT)
Number theoretic transform over a ring, also known as the discrete Fourier transform over a finite field, was introduced by Pollard [18]. For a finite field Fq and a generic sequence (a) of length d whose entries are from Fq , the forward DFT of (a) over Fq , denoted by (A), can be computed as Aj =
d−1 X
ai rij , 0 ≤ j ≤ d − 1 .
(1)
i=0
Here we refer to the elements of (a) and (A) by ai and Ai , respectively, for 0 ≤ i ≤ d − 1. Likewise, the inverse DFT of (A) over Fq can be computed as ai =
d−1 1 X · Aj r−ij , 0 ≤ i ≤ d − 1 . d
(2)
j=0
The above DFT computations over Fq are defined by utilizing a d-th primitive root of unity, denoted by r, from Fq or a finite extension of Fq . Definition 1 r is a primitive dth root of unity modulo n if rd = 1 (mod n) 3
and rd/t − 1 6= 0 (mod n) for any prime divisor t of d. Note that unlike the complex number r = ej2π/n used for the DFT computations in the complex number field, a finite field element r is used in finite fields. Also, the DFT over a finite field has no physical meaning as in the case of the DFT over the complex number field and serves only for the purpose of computing convolution efficiently as we will explain in section 2.3.2. However, in this paper we will still use the terms time domain and frequency domain for identifying a sequence and its number theoretic transform. Hence, for instance, we will refer to (a) and (A) as the time and frequency domain representations, respectively, of the same sequence. We would like to caution the reader that in a number theoretic transform the modulus q and the transform length d can not be chosen independently of each other. If the modulus is a composite number of the form q = pe11 pe22 pe33 . . . pekk , then the length of the transform, d, must divide gcd(p1 − 1, p2 − 1, p2 − 1, . . . , pk − 1). In the case of NTT over a finite field, i.e. when q is a prime p or a prime power pm , d must divide p − 1. 2.3.1
Mersenne Transform
Definition 2 A number of the form Mn = 2n −1, where n is a positive integer, is called a Mersenne number. If a Mersenne number is also a prime number, it is called a Mersenne prime. Mersenne primes are favourable choices as finite field characteristics due to their computational advantage in the modular reduction operation. Modular reduction by Mersenne primes can be achieved with simple addition and shift operations. A Number Theoretic Transform of special interest is Mersenne Transform, which is an NTT with arithmetic modulo a Mersenne number [19]. Mersenne Transform allows for very efficient forward and inverse transform operations for r = ±2. Multiplication of an n-bit number with integer powers of 2 modulo a Mersenne prime Mn can be performed by simple rotation of the bits of the n-bit number given as in the following properties. Property 1 For a positive integer i and an n-bit integer z ∈ GF (q), where q = 2n − 1 is a Mersenne prime, 2i · z mod q is equal to the n-bit number z rotated by i bits to the left. Property 2 For a positive integer i and an n-bit integer z ∈ GF (q), where q = 2n − 1 is a Mersenne prime, 2−i · z mod q is equal to the n-bit number z rotated by i bits to the right. Similarly, multiplication of an n-bit number with integer powers of −2 modulo a Mersenne prime Mn can be performed by simple rotation of the bits of the n-bit number, in addition to a negation if the power of −2 is odd. Hence, in a Mersenne transform with r = ±2 all of the multipllications by powers of r in the forward and inverse DFT computations can be achieved with simple bitwise rotations. In this case, for a transform length of d, the forward DFT computation can be achieved by only (d − 1)2 simple rotations and d(d − 1) additions/subtractions avoiding any multiplications. Likewise, the inverse DFT computation can be achieved by only (d − 1)2 simple rotations and d(d − 1) additions/subtractions, in addition to d constant multiplications by 1/d. For a more detailed complexity analysis of Mersenne transform, we refer the interested reader to [19]. Note that in a Mersenne transform the values of the sequence length d and the dth primitive root of unity r are dependent on each other. For the Mersenne transform with modulus q = Mn = 2n −1 and for r = ±2, the following equalities hold determining the relationship between d and r: 4
d=
n,
r=2.
2n ,
r = −2 .
For more information on several properties of Mersenne transform, the reader is referred to [19]. 2.3.2
Convolution Theorem and Multiplication in GF (q m ) Using NTT
A significant application of the Fourier transform is convolution. According to the convolution theorem, computing the convolution of two sequences in the time domain is equivalent to computing the componentwise parallel multiplication of their frequency domain representations. Convolution of any two d-element sequences (a) and (b) in the time domain results in another d-element sequence (c) and can be computed as follows: ci =
d−1 X
aj bi−j mod d , 0 ≤ i ≤ d − 1 .
(3)
j=0
The above convolution operation in the time domain is equivalent to the following computation in the frequency domain: Ci = Ai · Bi , 0 ≤ i ≤ d − 1 , (4) where (A), (B) and (C) denote the discrete Fourier transforms of (a), (b) and (c), respectively. Hence, convolution of two d-element sequences in the time domain, with complexity O(d2 ), is equivalent to simple pairwise multiplication of the DFTs of these sequences and has a surprisingly low O(d) complexity. Note in (3) that, the convolution is computed by making the sequences (a) and (b) periodic with d. This is equivalent to the cyclic convolution of the two sequences. Multiplication of two polynomials is basically the same as the acyclic (linear) convolution of the polynomial coefficients. We have seen that cyclic convolution can be performed very efficiently in the Fourier domain by pairwise coefficient multiplications. Hence, it will be wise to represent an element of GF (q m ), in polynomial representation an (m − 1)st degree polynomial with coefficients in GF (q), with at least a d = (2m − 1) element sequence by appending zeros at the end, so that the cyclic convolution of two such sequences will be equivalent to their acyclic convolution and give us their polynomial multiplication. We can form sequences by taking the ordered coefficients of polynomials. For instance, a(x) = a0 + a1 x + a2 x2 + . . . + am−1 xm−1 , an element of GF (q m ) in polynomial representation, can be interpreted as the following sequence after appending d − m zeros to the right: (a) = (a0 , a1 , a2 , . . . , am−1 , 0, 0, . . . , 0) .
(5)
For a(x), b(x) ∈ GF (q m ), and for d ≥ 2m − 1, the cyclic convolution of (a) and (b) yields a sequence (c) whose entries can be interpreted as the coefficients of a polynomial c(x) such that c(x) = a(x) · b(x). The computation of this cyclic convolution can be performed by simple pairwise coefficient multiplications in the discrete Fourier domain. The following straightforward algorithm realizes the polynomial multiplication c(x) = a(x) · b(x), without modular reduction by the field 5
generating polynomial, using the direct application of DFT. Algorithm for Polynomial Multiplication by Direct Application of DFT Input: a(x), b(x) ∈ GF (q m ) Output: c(x) = a(x) · b(x) Step 1. Interpret the coefficients of a(x) and b(x) as elements of the sequences (a) and (b), respectively, and append zeros to the right to make their lengths exactly d. Step 2. Convert (a) and (b) into their respective frequency domain representations (A) and (B) using the DFT operation of (1). Step 3. Multiply (A) and (B) together to compute (C) as in (4). Step 4. Convert (C) back to the time domain representation (c) using the inverse DFT operation of (2). Step 5. Interpret the first 2m − 1 coefficients of (c) as the coefficients of the product polynomial c(x). 2 The above method can be applied to multiplication of two polynomials of arbitrary degrees by only changing the condition d ≥ 2m − 1 to d > deg(a(x)) + deg(b(x)), where deg(p(x)) denotes the degree of p(x). Using this approach, polynomial multiplication is achieved with d(d − 1) additions and (d − 1)2 multiplications for the forward transform operation, d(d − 1) additions, (d − 1)2 multiplications and d constant multiplications for the inverse transform operation, and d multiplications for the multiplication in the frequency domain, totalling 2(d−1)2 +d multiplications, d constant multiplications and 2d(d − 1) additions in terms of coefficient arithmetic. However, as we mentioned in section 2.3.1, the same operations can be achieved with only 2(d − 1)2 bitwise rotations, d multiplications, d constant multiplications and 2d(d − 1) additions of coefficients when we use Mersenne transform. Note that, with this method the polynomial product c(x) = a(x) · b(x) is computed in the frequency domain but the final reduction by the field generating polynomial is not performed. One needs to convert (C) back to the time domain so that it can be reduced modulo the field generating polynomial and further multiplications can be performed on it using the same method. In the next section, we introduce a new method that performs both polynomial multiplication and modular reduction in the frequency domain and hence avoids costly back and forth conversions between the time and frequency domains.
3
Modular Multiplication in the Frequency Domain
In many finite field applications, a chain of arithmetic operations needs to be performed, rather than a solitary one. For example, in elliptic curve cryptography a scalar point product is computed by applying a chain of finite field additions, subtractions, multiplications, squarings and inversions on the input point coordinates [8]. The Montgomery residue representation has proven to be useful in this computation [16], [13], [12]. In using this method, first the operands are converted to their respective Montgomery residue representations, then utilizing Montgomery arithmetic the desired computation is implemented, and finally the result is converted back to the normal integer or polynomial representation. If there is a large number of operations performed in the Montgomery domain, due to the efficiency of the intervening computations, the forward and backward conversion operations become affordable. We introduce the same notion for the frequency domain. We want
6
to find an arithmetic operation in the frequency domain that is equivalent to ordinary modular or Montgomery multiplication in the time domain. For elliptic curve cryptography, we assume projective coordinates are used and all inversions are also performed by means of multiplications in the frequency domain. Thus, all finite field arithmetic operations, including modular reductions, may be performed in the frequency domain, avoiding any cost due to conversions between the time and the frequency domains. In the remainder of this section, we introduce the DFT modular multiplication algorithm which allows for both polynomial multiplication and modular reduction operations in the frequency domain. In section 3.1 we briefly present the notation we use in representing finite field elements in the frequency domain in the DFT modular multiplication algorithm and we present the algorithm in section 3.2.
3.1
Representation of Finite Field Elements in the DFT Modular Multiplication Algorithm
Let GF (q m ) denote a finite field extension generated by an irreducible polynomial f (x). For a(x), b(x) ∈ GF (q m ), both (m − 1)st degree polynomials with coefficients in GF (q) in polynomial representation, we let (a) and (b), respectively, denote their corresponding d-element time domain sequences. Note that the sequences (a) and (b) are extended to have length d by setting their higher d − m elements to 0. In this case, a(x) and b(x) have a maximum number of m coefficients each and, as explained in Sections 2 and 3, for the DFT operations we set d ≥ 2m − 1 and use r which is a d-th primitive root of unity in F = GF (q) or a finite extension of F. We denote the discrete Fourier transforms of (a) and (b) by (A) and (B), respectively, and use Ai and Bi , for 0 ≤ i ≤ d − 1, to denote their elements. Likewise, we represent the polynomial product c(x) = a(x) · b(x) with the d-element sequence (c). We let (C) represent the discrete Fourier transform of (c) and denote its elements by Ci , for 0 ≤ i ≤ d − 1. Finally, we let (f ) and (x) represent the d-element sequences for the irreducible field generating polynomial f (x) and the constant polynomial x ∈ GF (q m ), respectively, in the time domain. We denote the discrete Fourier transforms of (f ) and (x) by (F ) and (X), respectively, and let Fi and Xi , for 0 ≤ i ≤ d − 1, denote their elements. Hence, for instance, the constant polynomial x is represented by a d-element sequence as (x) = (0, 1, 0, 0, · · · , 0) , and the DFT of (x) is represented by the d-element sequence (X) given as (X) = (1, r, r2 , r3 , r4 , r5 , . . . , rd−1 ) . Note that we denote the elements of a sequence with the name of the sequence and a subscript for showing the location of the particular element in the sequence, e.g., X0 = 1 and Xd−1 = rd−1 .
3.2
The DFT Modular Multiplication Algorithm
In Table 1 we present the DFT modular multiplication algorithm which performs modular multiplication of a(x) and b(x), both elements of GF (q m ), in the discrete Fourier domain. We also give the complexity of each step of the algorithm in terms of the number of GF (q) multiplications and additions. The algorithm consists of two parts; multiplication (steps 1 and 2) and Montgomery reduction (steps 3 through 9). Multiplication is performed simply by pairwise multiplication of the elements of the frequency domain sequence representations of the two operands, i.e., Ai and Bi for 0 ≤ i ≤ d − 1. Multiplication of the two polynomials a(x) and b(x), each with degree at most (m − 1), results in another polynomial c(x) of degree at most 2m − 2. For performing further 7
multiplications over c(x) in the frequency domain, one needs to first reduce it modulo f (x) so that its time domain representation is of degree less than or equal to m − 1. Reduction is more complex and performed by Montgomery reduction in the frequency domain. The following algorithm achieves Montgomery reduction in our DFT modular multiplication algorithm. In the reduction process we use the normalized field generating polynomial f 0 (x) which equals f (x)/f0 mod q. Here, we let (f 0 ) denote the d-element time domain sequence representation for f 0 (x), with the sequence elements fi0 = fi · (f0−1 mod q) for 0 ≤ i ≤ d − 1. Hence, f 0 (x) is equivalent to f (x) but normalized to have its constant coefficient f 0 (0), or equivalently f00 in its time domain sequence representation (f 0 ), set to 1. Likewise, (F 0 ) denotes the DFT of (f 0 ) and its sequence elements are Fi0 = Fi · (f0−1 mod q), for 0 ≤ i ≤ d − 1. Note that, the reduction step of our algorithm results in the product c(x) of degree at most m − 1 by reducing the initial product modulo f 0 (x) which is equivalent to reducing it modulo f (x). However, it also adds an x−(m−1) factor to the otherwise correct result. DFT Montgomery Reduction Input: c(x) = a(x) · b(x) Output: c(x) = a(x) · b(x) · x−(m−1) mod f (x) Step 1. Compute −c0 , the negative of the first coefficient in the time domain sequence representation (c) of c(x), and put it in S (steps 4, 5, 6 and 7 in Table 1). Step 2. Make c0 zero by adding S · (f 0 ) to (c). Then shift the elements of new (c), which now has c0 = 0 and is a multiple of (x), by one position to the left via division by (x) (steps 8 and 9 in Table 1). Step 3. Repeat 1 and 2 for m − 1 times. Proof of Correctness: DFT Montgomery reduction is a direct adaptation of Montgomery reduction. We compute S such that (c(x) + S · f 0 (x)) is a multiple of x. We then divide (c(x) + S · f 0 (x)) by x and obtain a result which is congruent to c(x) · x−1 modulo f (x). By repeating this m − 1 times we obtain the final result which is congruent to c(x) · x−(m−1) modulo f (x). 2 The inputs of the DFT modular multiplication algorithm in Table 1 are the DFTs (A) and (B) of the d-element sequences (a) and (b) which represent the polynomials a(x) and b(x), respectively. The output of the algorithm is (C), the DFT of the sequence (c), where (c) represents the delement sequence for the polynomial c(x) = a(x)·b(x)·x−(m−1) mod f (x). The extra x−(m−1) factor shows that the DFT modular multiplication algorithm actually computes Montgomery products of polynomials [12] in the frequency domain. Hence, the input polynomials a(x) and b(x) may be viewed as the Montgomery residue representations of two polynomials u(x) and v(x) such that a(x) = u(x) · xm−1 mod f (x) and b(x) = v(x) · xm−1 mod f (x) . By the DFT modular multiplication the residue representation is kept intact, i.e.,
8
Step 1 Step 2 Step Step Step Step Step Step Step
3 4 5 6 7 8 9
MULTIPLICATION for i = 0 to d − 1 do Ci = Ai · Bi REDUCTION for j = 0 to m − 2 do S=0 for i = 0 to d − 1 do S = S + Ci S = −S/d for i = 0 to d − 1 do Ci = (Ci + Fi0 · S) · Xi−1
#Multiplication – d
#Add./Subtract. – –
– – – – m−1 – 2 · (m − 1) · d
– – – (m − 1) · d – – (m − 1) · d
Table 1: DFT modular multiplication algorithm for modular multiplication in GF (q m ) in the discrete Fourier domain
a(x) · b(x) · x−(m−1) mod f (x) = (u(x) · v(x)) · xm−1 mod f (x) which allows for further computation in the frequency domain. At first glance the complexity of the DFT modular multiplication algorithm appears to be even worse than that of the schoolbook multiplication method. For d ≈ 2m, modular multiplication in GF (q m ) with the DFT modular multiplication method requires 4m2 − m − 1 multiplications and 4m2 − 4m additions in the ground field GF (q), while the schoolbook method requires only m2 multiplications and (m − 1)2 additions ignoring the cost of modular reduction. Nevertheless, this complexity may be improved dramatically by using special values for q, r, d and the irreducible field generating polynomial f (x) as we will explain in the next section.
4
Efficient Implementation of DFT Modular Multiplication over Finite Fields Constructed with Special Parameters
In this section, we investigate the existence of finite fields of the form GF (q m ), where q is a Mersenne prime, and efficient parameters such as irreducible field generating polynomials of special form, sequence lengths satisfying d ≥ 2m − 1 and d ≈ 2m, and d-th primitive roots of unity r of special form exist for efficient implementation of DFT modular multiplication. We would like to caution the reader that all these parameters are dependent on each other and can not be chosen independently.
4.1
Sequence Length d
In section 2.3.2 we have seen that a necessary condition for performing polynomial multiplication in the frequency domain is that d ≥ 2m − 1. Using the smallest possible sequence length d will lead to smallest number of arithmetic operations in the computation of DFT modular multiplication (see Table 1). Optimally, d = 2m − 1 will lead to the least number of arithmetic operations. Furthermore, note that multiplications by −d−1 in step 7 of the DFT modular multiplication algorithm are constant multiplications, and hence can be computed by a sequence of additions and shift operations. For an n-bit field characteristic q, these constant multiplications can be achieved with n/2 shifts and n/2 additions, on average. And naturally, for lower Hamming weight 9
−d−1 values, the constant multiplications by −d−1 can be performed more efficiently with lower number of shifts and additions. This would be very helpful especially for implementations on processors without a multiplication instruction where multiplications are performed with shift and add instructions. However, again we would like to note that the parameter d can not be selected independently and its value depends on other parameters such as the field characteristic q and the dth primitive root of unity r.
4.2
Mersenne Primes as Finite Field Characteristics
As mentioned in section 2.3.1, the use of a Mersenne number as the modulus and selection of r = ±2 allow for very efficient DFT computations in Mersenne transform. We will utilize this fact in our DFT modular multiplication algorithm. In a finite field with a Mersenne prime characteristic, a multiplication by ri , for r = ±2 and i an integer, can be achieved by a simple bitwise rotation (in addition to a negation when r = −2 and i is odd) which is inexpensive (see section 2.3.1). In our DFT modular multiplication algorithm of Table 1, this property may be exploited if we choose q to be a Mersenne prime. In step 9 of the algorithm, Xi−1 = r−i and hence multiplications by Xi−1 can be computed with simple bitwise rotations. Multiplications by Fi0 can also be avoided in a similar fashion. For instance, for f (x) = xm + rs0 with s0 an integer, Fi0 = rmi−s0 mod d + 1 ,
(6)
and hence for r = ±2 multiplications by Fi0 can be performed with only one bitwise rotation and one addition/subtraction.
4.3
Binomials or Trinomials as Field Generating Polynomials
The efficiency of step 9 of the DFT modular multiplication algorithm can be improved significantly by using efficient parameters for r and the field generating polynomial f (x). The critical operations in step 9 are the two multiplications, which are multiplications of Fi0 by S and Ci + Fi0 · S by Xi−1 . The multipliers Fi0 are constants and their values depend on the irreducible polynomial f (x), as well as on the value of r. Similarly, the multipliers Xi−1 are constants whose values depend on the value of r. When r = 2, Xi−1 becomes r−i = 2−i = 2d−i and hence when the finite field characteristic is a Mersenne prime multiplications by Xi−1 , as in step 9 of the DFT modular multiplication algorithm, become simple (d − i)-bit left-rotations (or equivalently i-bit right-rotations), which have negligible cost compared to regular multiplication. Similarly, when r = −2, multiplications by Xi−1 = r−i = (−2)−i = (−2)d−i become simple (d − i)-bit left-rotations (or equivalently i-bit right-rotations) and a negation depending on whether the power of −2 is odd. The multiplications by Fi0 in step 9 of the DFT modular multiplication algorithm can be simplified significantly by using a low Hamming weight field generating polynomial for f (x). We can guarantee that the elements of the sequence (F ) are in the form of summations of powers of r by selecting the coefficients of f (x) to be 1 or, for that matter, any integer power of r. When f0 and all other coefficients of f (x) are powers of r, so are the coefficients of f 0 (x) = f (x)/f0 . In this case, the elements of the sequence (F 0 ) will also be in the form of summations of powers of r. Hence, for r = ±2 and in a finite field with a Mersenne prime characteristic, multiplications by Fi0 can be performed with simple bitwise rotations and additions/subtractions. For clarification, we make the following illustration. For an irreducible field generating polynomial of the form f (x) = xm + rsm−1 xm−1 + rsm−2 xm−2 + rsm−3 xm−3 + · · · + rs2 x2 + rs1 x + rs0 , 10
whose coefficients are all integer powers of r, the d-element time domain sequence representation is (f ) = (rs0 , rs1 , rs2 , rs3 , · · · , rsm−2 , rsm−1 , 1, 0, 0, · · · , 0) and the DFT of (f ) is (F ) = (F0 , F1 , F2 , F3 , · · · , Fd−1 ) , where Fi = rmi + r(m−1)i+sm−1 + r(m−2)i+sm−2 + r(m−3)i+sm−3 + · · · + r2i+s2 + ri+s1 + rs0 . Likewise, for the normalized field generating polynomial f 0 (x) = f (x)/f0 = r−s0 xm + rsm−1 −s0 xm−1 + rsm−2 −s0 xm−2 + rsm−3 −s0 xm−3 + · · · + rs1 −s0 x + 1 , the d-element time domain sequence representation is (f 0 ) = (1, rs1 −s0 , rs2 −s0 , rs3 −s0 , · · · , rsm−2 −s0 , rsm−1 −s0 , r−s0 , 0, 0, · · · , 0) and the DFT of (f 0 ) is 0 (F 0 ) = (F00 , F10 , F20 , F30 , · · · , Fd−1 ),
where Fi0 = Fi /f0 = rmi−s0 + r(m−1)i+sm−1 −s0 + r(m−2)i+sm−2 −s0 + r(m−3)i+sm−3 −s0 + · · · + ri+s1 −s0 + 1 . Note that all elements of the sequence (F 0 ), namely Fi0 for 0 ≤ i ≤ d − 1, are in the form of summations of powers of r. Therefore, when a finite field with a Mersenne prime characteristic is used and for r = ±2, multiplications by Fi0 can be performed with simple bitwise rotations and additions/subtractions. Clearly, selecting low Hamming weight polynomials for f (x) may reduce the complexity of multiplications by Fi0 by means of reducing the number of required additions/subtractions, as we shall see better next. In Table 2 we list the DFTs for some simple polynomials. Our intention is to determine the type of irreducible polynomial that has a simpler DFT and therefore works best for our purpose. Since the DFT is a linear operation, the DFT of the sum of two polynomials equals the sum of their individual DFTs. Similarly, the DFT of a constant multiple of a polynomial equals the same constant multiple of the DFT of that polynomial. Hence, the DFT of any polynomial is a summation of constant multiples of sequences such as those given in Table 2. f (x) r x xm
(F ) ( r ( 1 ( 1
, , ,
r r rm
, , ,
r r2 r2m
, , ,
r r3 r3m
, , ,
r r4 r4m
, , ,
... ... ...
, , ,
r ) rd−1 ) r(d−1)m )
Table 2: Frequency domain sequence representations for some simple polynomials Naturally, an irreducible binomial of the form f (x) = xm ± rs0 with an integer s0 , or trinomial 0 of the form f (x) = xm ± rsm0 xm ± rs0 with m0 < m and integers sm0 and s0 , will have a simpler frequency domain sequence representation and multiplications by F 0 , as in step 9 of the DFT 11
modular multiplication algorithm in Table 1, will be easier to perform. For instance, provided that f (x) = xm + rs0 we will have Fi = rmi + rs0 and Fi0 = f0−1 Fi = rmi−s0 + 1 . And for r = ±2 and a Mersenne prime q, multiplication of a number in GF (q) by Fi0 , of the form rmi−s0 + 1, can easily be performed with only a bitwise rotation and an addition/subtraction. 0 Likewise, for the trinomial f (x) = xm + rsm0 xm + rs0 , 0
Fi0 = rmi−s0 + rm i+sm0 −s0 + 1 , and multiplications by Fi0 can be performed with only two bitwise rotations and two additions/subtractions. Also, note that for any integer s the equality rs = rs mod d holds, since r is a dth primitive root of unity in GF (q), a condition satisfied for the application of the DFT, and therefore all required bitwise rotations due to multiplications by powers of r will be for less than d bits. Moreover, since r is a dth primitive root of unity, all numbers of the form rs , for any integer s, will have only d distinct values. Hence, in step 9 of the DFT modular multiplication algorithm, all but d of the multiplications of S by rk , due to multiplications by Fi0 and Xi−1 , may be avoided by precomputing the d distinct values of the rk · S products, for the value of S computed at step 7, and using these precomputed values as needed in the computations of Ci = (Ci + Fi0 · S) · Xi−1 .
4.4
Existence of Efficient Parameters
In Table 3 we give a list of parameters such as Mersenne prime Mn , prime extension degree m, sequence length d, r = ±2 and equivalent binary field size for which efficient modular multiplication in GF (q m ) using the DFT modular multiplication algorithm is possible for elliptic curve cryptography. Note that between d and n we always have the relationships d = n or d = 2n for r = 2 or r = −2, respectively. Also note that in our list we include only those cases where the field extension degree m is prime, which is a security requirement for elliptic curve cryptography over extension degrees. For each case listed in Table 3, in Appendix we also provide a list of irreducible binomials or trinomials as field generating polynomials that would allow for efficient implementation of DFT modular multiplication. For some of the cases we were able to find efficient irreducible binomials and included them, while for other cases we were not able to find such binomials and included trinomials instead. However, we would like to note that for the computationally more desirable cases of m = d/2 we were always able to find efficient irreducible binomials.
5
Complexity of the DFT Modular Multiplication Algorithm
In this section, we present the complexity of the DFT modular multiplication algorithm for a practical set of parameters and compare it with three other efficient multiplication approaches, namely the direct application of DFT multiplication with time domain modular reduction (see section 2.3.2), algorithm by Bajard et al. [3] and Karatsuba algorithm [11]. In our complexity 12
n 13 17 17 17 19 19 19 19 31 31
q = Mn 8191 131071 131071 131071 524287 524287 524287 524287 2147483647 2147483647
m 13 11 13 17 11 13 17 19 11 13
d 26 34 34 34 38 38 38 38 31 31
r −2 −2 −2 −2 −2 −2 −2 −2 2 2
equivalent binary field size ∼ 2169 ∼ 2187 ∼ 2221 ∼ 2289 ∼ 2209 ∼ 2247 ∼ 2323 ∼ 2361 ∼ 2341 ∼ 2403
Table 3: List of Mn , m, d and r = ±2 values that would make efficient DFT modular multiplication in GF (q m ) possible for elliptic curve cryptography over finite fields of size 150 to 500 bits.
1 2 3 4 5 6 7 8 9
MULTIPLICATION for i = 0 to d − 1 do Ci = Ai · Bi REDUCTION for j = 0 to m − 2 do S=0 for i = 0 to d − 1 do S = S + Ci S = −S/d for i = 0 to d − 1 do Ci = (Ci + Fi0 · S) · Xi−1
#Mult. – d ≈ 2m
#Const. Mult. – –
#Add./Subtract. – –
#Rotation – –
– – – –
–
– – – –
– – –
m−1 – –
– – – (m − 1) · d ≈ 2m(m − 1) – – 2(m − 1) · d ≈ 4m(m − 1)
– –
– – 2 · (m − 1) · d ≈ 4m(m − 1)
Table 4: DFT modular multiplication algorithm in GF (q m ) and its complexity for a Mersenne prime q, an irreducible field generating binomial of the form f (x) = xm + rs0 with an integer s0 , r = ±2 and d ≈ 2m.
analysis we assume the use of efficient parameters such as Mersenne prime finite field characteristics as q, irreducible field generating binomials of the form f (x) = xm + 2s0 with a positive integer s0 , d-th primitive root of unity r = ±2 and sequence length d ≈ 2m. The field parameters we use, such as low Hamming weight field generating polynomial and Mersenne prime field characteristics, would lead to efficient implementation of multiplication for all four methods. Therefore, for the selected parameters, we can safely assume that our comparisons are fair. In our complexity derivations for the DFT modular multiplication we utilize all efficient implementation ideas mentioned in Section 4 and obtain the complexity of multiplication in GF (q m ) in terms of the number of GF (q) operations as shown in Table 4. Note that, for performing DFT modular multiplication in GF (q m ) only a linear number of GF (q) multiplications (step 2) and constant multiplications (step 7) are needed, in addition to some additions/subtractions (steps 6, 9) and rotation operations (step 9) in GF (q) which are relatively inexpensive. In Table 5, we present the complexities of GF (q m ) multiplication in terms of GF (q) operations for all four approaches. The complexity presented for Karatsuba algorithm is only an approximate one for the case when m is a power of two. For the best complexities of Karatsuba algorithm for 13
#Multiplication #Const. Mult. #Add./Subtr. #Rotation
DFT Multiplication (direct approach) d ≈ 2m d ≈ 2m 2d(d − 1) ≈ 8m2 − 4m 2(d − 1)2 ≈ 8m2 − 8m + 2
Bajard et al.
Karatsuba
3m 2m + 2m
≈ mlog2 3 –
2m2 − m
≈ 6mlog2 3 − 7m + 1
–
m−1
2
DFT Modular Multiplication d ≈ 2m m−1 3(m − 1)d ≈ 6m2 − 6m 2(m − 1)d ≈ 4m2 − 4m
Table 5: Complexities of modular multiplication in GF (q m ) where the field generating polynomial is f (x) = xm + 2s0 for an integer s0 , q is a Mersenne prime and d ≈ 2m
different extension degrees, the reader is referred to [22]. Note that the number of ground field multiplications required by Karatsuba algorithm is subquadratic, whereas all other approaches require only linear number of ground field multiplications. Furthermore, although the algorithm by Bajard et al. has only linear complexity in terms of the number of multiplications, it still requires performing a quadratic number of constant multiplications in GF (q). Whereas the number of constant multiplications required by DFT modular multiplication is only linear. In the number of additions and subtractions, Karatsuba algorithm has subquadratic complexity whereas the complexities of the other three methods are quadratic. Finally, multiplication with the direct DFT approach and DFT modular multiplication algorithm require a quadratic number of simple bitwise rotations, whereas Karatsuba multiplication requires only a linear number of bitwise rotations and the method by Bajard et al. does not require any. Clearly, the complexity of the DFT modular multiplication algorithm is better than the direct DFT approach. Moreover, since DFT modular multiplication requires significantly less number of complex operations such as multiplication and constant multiplication, its overall performance appears to be better than the algorithm by Bajard et al. and the Karatsuba algorithm, especially for computational platforms where multiplication is expensive compared to other operations such as addition, subtraction or bitwise rotation. In many modern microprocessors speed of a multiplication instruction is comparable to or only a few times slower than an addition instruction. However this is not the case especially for some microprocessors where there are no special multiplication instructions and multiplications are performed with a sequence of shifts and additions. For instance, ATMEL 90LS8535 [1] is a microprocessor that is widely used in embedded applications, yet it has no multiplication instructions. Likewise, in hardware complexity of multiplication is significantly more than complexity of addition. For instance, a straightforward implementation of an n-bit multiplication can be achieved via n additions and n shift operations. Therefore, in serialized implementation the complexity of an n-bit multiplication may be assumed to be roughly equal to the complexity of n n-bit additions and n shift operations. Using the same approach, we may assume that multiplication by a constant n-bit number can be performed with n/2 shifts and n/2 additions on average, and hence has complexity equal to n/2 addition and n/2 shift operations. Under these assumptions, Table 5 gives the complexities of modular multiplication in GF (q 13 ) for all four methods where q = 213 − 1 and GF (q 13 ) is constructed using the irreducible binomial f (x) = x13 − 2. The table also includes the total number of clock cycles that a single multiplication with these methods takes, assuming addition/subtraction and rotation operations all take a single clock cycle. Note that this finite field has size ∼ 2169 and can be considered suitable for elliptic curve cryptography. We caution the reader that all values in Table 5 are obtained directly from
14
#Multiplication #Const. Mult. #Add./Subtr. #Rotation #Total Clock Cycles
DFT Multiplication (direct approach) 26 26 1300 1250 3564
Bajard et al.
Karatsuba
39 364 325 – 6071
91 – 390 12 2768
DFT Modular Multiplication 26 12 936 624 2392
Table 6: Complexities of modular multiplication in GF (q 13 ) where the field generating polynomial is f (x) = x13 − 2, q = 213 − 1, d = 26 and r = −2
Table 5, except the complexity of Karatsuba multiplication for GF (q 13 ) which we draw from [22] . Also, note that for all four methods we assume the average case for the complexities of constant multiplications regardless of the values of the particular parameters used. Hence, as demonstrated with the examples above, finite field polynomial multiplication may be achieved more efficiently with DFT modular multiplication than with other efficient methods such as the Karatsuba algorithm in computational platforms where multiplication operation is significantly more expensive than simple operations such as addition, subtraction and bitwise rotation.
6
Conclusion
We introduced the DFT modular multiplication algorithm which performs modular multiplication in the frequency domain using Montgomery reduction. The FFT based multiplication technique, which performs multiplication in the frequency domain, is known to be very efficient. Nevertheless, due to the lack of a frequency domain reduction algorithm, for performing modular reductions conversion to the time domain is needed, which adds a significant overhead and makes the FFT based multiplication algorithms impractical for small operands. In this work, by allowing for modular reductions in the frequency domain, our algorithm avoids back and forth conversions between the frequency and the time domains. We have shown that, with our method, in computational platforms where multiplication is expensive frequency domain modular multiplication for finite fields becomes practical and can be achieved more efficiently than multiplication in the time domain for finite field sizes relevant to elliptic curve cryptography.
References [1] Atmel AVR 8bit RISC processor. http://www.atmel.com/dyn/resources/prod documents/1041S.PDF. [2] J.-C. Bajard and L. Imbert. A Full RNS Implementation of RSA. IEEE Transactions on Computers, 53(6):769–774, June 2004. [3] J.-C. Bajard, L. Imbert, C. N`egre, and T. Plantard. Efficient Multiplication in GF (pk ) for Elliptic Curve Cryptography. In Proceedings of the 16th IEEE Symposium on Computer Arithmetic (ARITH’03). IEEE, 2003. [4] E. R. Berlekamp. Algebraic Coding Theory. McGraw-Hill, New York, New York, USA, 1968. 15
[5] R. E. Blahut. Theory and Practice of Error Control Codes. Addison-Wesley, Reading, Massachusetts, USA, 1983. [6] R. E. Blahut. Fast Algorithms for Digital Signal Processing. Addison-Wesley, Reading, Massachusetts, USA, 1985. [7] I.F. Blake, X.H. Gao, R.C. Mullin, S.A. Vanstone, and T. Yaghgoobin. Applications of Finite Fields. Kluwer Academic, 1999. [8] I.F. Blake, G. Seroussi, and N. Smart. Elliptic Curves in Cryptography. Cambridge University Press, London Mathematical Society Lecture Notes Series 265, 1999. [9] J. Cooley and J. Tukey. An Algorithm for the Machine Calculation of Complex Fourier Series. Mathematics of Computation, 19:297–301, 1965. [10] R. Crandall and C. Pomerance. Prime Numbers. Springer-Verlag, New York, NY, USA, 2001. [11] A. Karatsuba and Y. Ofman. Multiplication of Multidigit Numbers on Automata. Sov. Phys. Dokl. (English translation), 7(7):595–596, 1963. [12] C ¸ . K Ko¸c and T. Acar. Montgomery Multplication in GF (2k ). Design, Codes, and Cryptography, 14(1):57–69, 1998. ¸ . K Ko¸c, T. Acar, and B. Kaliski. Analyzing and Comparing Montgomery Multiplication [13] C Algorithms. IEEE Micro, pages 26–33, June 1996. [14] R. Lidl and H. Niederreiter. Finite Fields, volume 20 of Encyclopedia of Mathematics and its Applications. Addison-Wesley, Reading, Massachusetts, USA, 1983. [15] R. J. McEliece. Finite Fields for Computer Scientists and Engineers. Kluwer Academic Publishers, 2nd edition, 1989. [16] P. L. Montgomery. Modular Multiplication without Trial Division. Mathematics of Computation, 44(170):519–521, April 1985. [17] C. Paar. Efficient VLSI Architectures for Bit-Parallel Computation in Galois Fields. PhD thesis, (Engl. transl.), Institute for Experimental Mathematics, University of Essen, Essen, Germany, June 1994. ISBN 3–18–332810–0. [18] J. M. Pollard. The Fast Fourier Transform in a Finite Field. Mathematics of Computation, 25:365–374, 1971. [19] C. M. Rader. Discrete Convolutions via Mersenne Transforms. IEEE Transactions on Computers, C-21(12):1269–1273, December 1972. [20] R. L. Rivest, A. Shamir, and L. Adleman. A Method for Obtaining Digital Signatures and Public-Key Cryptosystems. Communications of the ACM, 21(2):120–126, February 1978. [21] A. Sch¨ onhage and V. Strassen. Schnelle Multiplikation großer Zahlen. Computing, 7:281–292, 1971. [22] A. Weimerskirch and C. Paar. Generalizations of the Karatsuba Algorithm for Efficient Implementations. Technical report, Department of Electrical Engineering and Information Sciences, Ruhr-Universit¨at Bochum, Germany. 16
Appendix x13 − 2 x13 − 22 x13 − 23 x13 − 24 x13 − 25 x13 − 26
x13 − 27 x13 − 28 x13 − 29 x13 − 210 x13 − 211 x13 − 212
x13 + 2 x13 + 22 x13 + 23 x13 + 24 x13 + 25 x13 + 26
x13 + 27 x13 + 28 x13 + 29 x13 + 210 x13 + 211 x13 + 212
Table 7: Irreducible binomials of the form x13 ±2s0 for 0 ≤ s0 ≤ 12 that can be used for construction of the finite field GF (q 13 ) where q = 213 − 1.
x11 + 28 x3 + 1 x11 + 211 x3 + 22 x11 + 22 x3 + 24 x11 + x3 + 26 x11 + 23 x3 + 28 x11 + 26 x3 + 210 x11 + 29 x3 + 212 x11 + x3 + 214 x11 + 23 x3 + 216
x11 + 213 x3 + 1 x11 + 216 x3 + 22 x11 + 214 x3 + 24 x11 + 25 x3 + 26 x11 + 28 x3 + 28 x11 + 211 x3 + 210 x11 + 214 x3 + 212 x11 + 212 x3 + 214 x11 + 215 x3 + 216
x11 + 2x3 + 2 x11 + 24 x3 + 23 x11 + 27 x3 + 25 x11 + 210 x3 + 27 x11 + 2x3 + 29 x11 + 24 x3 + 211 x11 + 22 x3 + 213 x11 + 25 x3 + 215
x11 + 26 x3 + 2 x11 + 29 x3 + 23 x11 + 212 x3 + 25 x11 + 215 x3 + 27 x11 + 213 x3 + 29 x11 + 216 x3 + 211 x11 + 27 x3 + 213 x11 + 210 x3 + 215
Table 8: Irreducible trinomials of the form x11 + 2s3 x3 + 2s0 for 0 ≤ s3 , s0 ≤ 16 that can be used for construction of the finite field GF (q 11 ) where q = 217 − 1.
17
x13 + 2x + 1 x13 + 212 x + 22 x13 + 26 x + 24 x13 + x + 26 x13 + 211 x + 28 x13 + 25 x + 210 x13 + x + 212 x13 + 210 x + 214 x13 + 24 x + 216
x13 + 22 x + 1 x13 + 213 x + 22 x13 + 27 x + 24 x13 + 2x + 26 x13 + 212 x + 28 x13 + 26 x + 210 x13 + 216 x + 212 x13 + 211 x + 214 x13 + 25 x + 216
x13 + 215 x + 2 x13 + 29 x + 23 x13 + 23 x + 25 x13 + 214 x + 27 x13 + 28 x + 29 x13 + 22 x + 211 x13 + 213 x + 213 x13 + 27 x + 215
x13 + 216 x + 2 x13 + 210 x + 23 x13 + 24 x + 25 x13 + 215 x + 27 x13 + 29 x + 29 x13 + 23 x + 211 x13 + 214 x + 213 x13 + 28 x + 215
Table 9: Irreducible trinomials of the form x13 + 2s1 x + 2s0 for 0 ≤ s1 , s0 ≤ 16 that can be used for construction of the finite field GF (q 13 ) where q = 217 − 1.
x17 − 2 x17 − 22 x17 − 23 x17 − 24 x17 − 25 x17 − 26 x17 − 27 x17 − 28
x17 − 29 x17 − 210 x17 − 211 x17 − 212 x17 − 213 x17 − 214 x17 − 215 x17 − 216
x17 + 2 x17 + 22 x17 + 23 x17 + 24 x17 + 25 x17 + 26 x17 + 27 x17 + 28
x17 + 29 x17 + 210 x17 + 211 x17 + 212 x17 + 213 x17 + 214 x17 + 215 x17 + 216
Table 10: Irreducible binomials of the form x17 ± 2s0 for 0 ≤ s0 ≤ 16 that can be used for construction of the finite field GF (q 17 ) where q = 217 − 1.
x11 + 29 x + 1 x11 + 2x + 22 x11 + 24 x + 24 x11 + 211 x + 26 x11 + 23 x + 28 x11 + 26 x + 210 x11 + 213 x + 212 x11 + 2x + 214 x11 + 28 x + 216 x11 + x + 218
x11 + 213 x + 1 x11 + 216 x + 22 x11 + 28 x + 24 x11 + 215 x + 26 x11 + 218 x + 28 x11 + 210 x + 210 x11 + 217 x + 212 x11 + 25 x + 214 x11 + 212 x + 216 x11 + 215 x + 218
x11 + 23 x + 2 x11 + 210 x + 23 x11 + 22 x + 25 x11 + 25 x + 27 x11 + 212 x + 29 x11 + x + 211 x11 + 27 x + 213 x11 + 214 x + 215 x11 + 22 x + 217
x11 + 27 x + 2 x11 + 214 x + 23 x11 + 217 x + 25 x11 + 29 x + 27 x11 + 216 x + 29 x11 + 24 x + 211 x11 + 211 x + 213 x11 + 218 x + 215 x11 + 26 x + 217
Table 11: Irreducible trinomials of the form x11 + 2s1 x + 2s0 for 0 ≤ s1 , s0 ≤ 18 that can be used for construction of the finite field GF (q 11 ) where q = 219 − 1.
18
x13 + 24 x + 1 x13 + 215 x + 24 x13 + 27 x + 28 x13 + 218 x + 212 x13 + 210 x + 216
x13 + 22 x + 2 x13 + 213 x + 25 x13 + 25 x + 29 x13 + 216 x + 213 x13 + 28 x + 217
x13 + x + 22 x13 + 211 x + 26 x13 + 23 x + 210 x13 + 214 x + 214 x13 + 26 x + 218
x13 + 217 x + 23 x13 + 29 x + 27 x13 + 2x + 211 x13 + 212 x + 215
Table 12: Irreducible trinomials of the form x13 + 2s1 x + 2s0 for 0 ≤ s1 , s0 ≤ 18 that can be used for construction of the finite field GF (q 13 ) where q = 219 − 1.
x17 + 25 x + 1 x17 + 211 x + 24 x17 + 217 x + 28 x17 + 24 x + 212 x17 + 210 x + 216
x17 + 216 x + 2 x17 + 23 x + 25 x17 + 29 x + 29 x17 + 215 x + 213 x17 + 22 x + 217
x17 + 28 x + 22 x17 + 214 x + 26 x17 + 2x + 210 x17 + 27 x + 214 x17 + 213 x + 218
x17 + x + 23 x17 + 26 x + 27 x17 + 212 x + 211 x17 + 218 x + 215
Table 13: Irreducible trinomials of the form x17 + 2s1 x + 2s0 for 0 ≤ s1 , s0 ≤ 18 that can be used for construction of the finite field GF (q 17 ) where q = 219 − 1.
x19 − 2 x19 − 22 x19 − 23 x19 − 24 x19 − 25 x19 − 26 x19 − 27 x19 − 28 x19 − 29
x19 − 210 x19 − 211 x19 − 212 x19 − 213 x19 − 214 x19 − 215 x19 − 216 x19 − 217 x19 − 218
x19 + 2 x19 + 22 x19 + 23 x19 + 24 x19 + 25 x19 + 26 x19 + 27 x19 + 28 x19 + 29
x19 + 210 x19 + 211 x19 + 212 x19 + 213 x19 + 214 x19 + 215 x19 + 216 x19 + 217 x19 + 218
Table 14: Irreducible binomials of the form x19 ± 2s0 for 0 ≤ s0 ≤ 18 that can be used for construction of the finite field GF (q 19 ) where q = 219 − 1.
19
x11 + 220 x + 1 x11 + 219 x + 22 x11 + 218 x + 24 x11 + 217 x + 26 x11 + 216 x + 28 x11 + 215 x + 210 x11 + 214 x + 212 x11 + 213 x + 214 x11 + 212 x + 216 x11 + 211 x + 218 x11 + 210 x + 220 x11 + 29 x + 222 x11 + 28 x + 224 x11 + 27 x + 226 x11 + 26 x + 228 x11 + 25 x + 230
x11 + 229 x + 1 x11 + 228 x + 22 x11 + 227 x + 24 x11 + 226 x + 26 x11 + 225 x + 28 x11 + 224 x + 210 x11 + 223 x + 212 x11 + 222 x + 214 x11 + 221 x + 216 x11 + 220 x + 218 x11 + 219 x + 220 x11 + 218 x + 222 x11 + 217 x + 224 x11 + 216 x + 226 x11 + 215 x + 228 x11 + 214 x + 230
x11 + 24 x + 2 x11 + 23 x + 23 x11 + 22 x + 25 x11 + 2x + 27 x11 + x + 29 x11 + 28 x + 211 x11 + 27 x + 213 x11 + 26 x + 215 x11 + 25 x + 217 x11 + 24 x + 219 x11 + 23 x + 221 x11 + 22 x + 223 x11 + 2x + 225 x11 + x + 227 x11 + 221 x + 229
x11 + 213 x + 2 x11 + 212 x + 23 x11 + 211 x + 25 x11 + 210 x + 27 x11 + 29 x + 29 x11 + 230 x + 211 x11 + 229 x + 213 x11 + 228 x + 215 x11 + 227 x + 217 x11 + 226 x + 219 x11 + 225 x + 221 x11 + 224 x + 223 x11 + 223 x + 225 x11 + 222 x + 227 x11 + 230 x + 229
Table 15: Irreducible trinomials of the form x11 + 2s1 x + 2s0 for 0 ≤ s1 , s0 ≤ 30 that can be used for construction of the finite field GF (q 11 ) where q = 231 − 1.
20
x13 + 28 x + 1 x13 + 2x + 2 x13 + 22 x + 22 x13 + 26 x + 23 x13 + 22 x + 24 x13 + x + 25 x13 + 24 x + 26 x13 + x + 27 x13 + 213 x + 28 x13 + 22 x + 29 x13 + 27 x + 210 x13 + 211 x + 211 x13 + x + 212 x13 + 25 x + 213 x13 + 29 x + 214 x13 + 22 x + 215 x13 + 23 x + 216 x13 + 27 x + 217 x13 + x + 218 x13 + 2x + 219 x13 + 25 x + 220 x13 + 2x + 221 x13 + 214 x + 222 x13 + 23 x + 223 x13 + 28 x + 224 x13 + 212 x + 225 x13 + 2x + 226 x13 + 26 x + 227 x13 + 210 x + 228 x13 + 23 x + 229 x13 + 24 x + 230
x13 + 212 x + 1 x13 + 24 x + 2 x13 + 217 x + 22 x13 + 210 x + 23 x13 + 211 x + 24 x13 + 215 x + 25 x13 + 28 x + 26 x13 + 29 x + 27 x13 + 217 x + 28 x13 + 26 x + 29 x13 + 222 x + 210 x13 + 215 x + 211 x13 + 24 x + 212 x13 + 220 x + 213 x13 + 213 x + 214 x13 + 25 x + 215 x13 + 218 x + 216 x13 + 211 x + 217 x13 + 23 x + 218 x13 + 216 x + 219 x13 + 29 x + 220 x13 + 210 x + 221 x13 + 218 x + 222 x13 + 27 x + 223 x13 + 223 x + 224 x13 + 216 x + 225 x13 + 25 x + 226 x13 + 221 x + 227 x13 + 214 x + 228 x13 + 26 x + 229 x13 + 219 x + 230
x13 + 215 x + 1 x13 + 213 x + 2 x13 + 221 x + 22 x13 + 213 x + 23 x13 + 226 x + 24 x13 + 219 x + 25 x13 + 211 x + 26 x13 + 224 x + 27 x13 + 220 x + 28 x13 + 29 x + 29 x13 + 226 x + 210 x13 + 218 x + 211 x13 + 27 x + 212 x13 + 224 x + 213 x13 + 216 x + 214 x13 + 214 x + 215 x13 + 222 x + 216 x13 + 214 x + 217 x13 + 212 x + 218 x13 + 220 x + 219 x13 + 212 x + 220 x13 + 225 x + 221 x13 + 221 x + 222 x13 + 210 x + 223 x13 + 227 x + 224 x13 + 219 x + 225 x13 + 28 x + 226 x13 + 225 x + 227 x13 + 217 x + 228 x13 + 215 x + 229 x13 + 223 x + 230
x13 + 224 x + 1 x13 + 228 x + 2 x13 + 224 x + 22 x13 + 222 x + 23 x13 + 230 x + 24 x13 + 222 x + 25 x13 + 220 x + 26 x13 + 228 x + 27 x13 + 229 x + 28 x13 + 218 x + 29 x13 + 229 x + 210 x13 + 227 x + 211 x13 + 216 x + 212 x13 + 227 x + 213 x13 + 225 x + 214 x13 + 229 x + 215 x13 + 225 x + 216 x13 + 223 x + 217 x13 + 227 x + 218 x13 + 223 x + 219 x13 + 221 x + 220 x13 + 229 x + 221 x13 + 230 x + 222 x13 + 219 x + 223 x13 + 230 x + 224 x13 + 228 x + 225 x13 + 217 x + 226 x13 + 228 x + 227 x13 + 226 x + 228 x13 + 230 x + 229 x13 + 226 x + 230
Table 16: Irreducible trinomials of the form x13 + 2s1 x + 2s0 for 0 ≤ s1 , s0 ≤ 30 that can be used for construction of the finite field GF (q 13 ) where q = 231 − 1.
21
1
Introduction
Finite fields have many applications in coding theory [5, 4], cryptography [15, 7, 17], and digital signal processing [6]. Hence efficient implementation of finite field arithmetic operations is crucial. Multiplication of field elements is commonly implemented in terms of modular multiplication of polynomials realized in two steps: ordinary polynomial multiplication and modular reduction of the result by the field generating polynomial. Unlike polynomial multiplication, modular reduction has only linear complexity, i.e. O(m), for polynomials of degree m − 1 when a fixed special modulus is chosen. Hence it is essential to attack the complexity of the multiplication step. The classical schoolbook polynomial multiplication method for multiplication of (m − 1)st degree polynomials has quadratic complexity, i.e. O(m2 ), given in terms of coefficient multiplications and additions. The complexity may be improved to O(mlog2 3 ) using the Karatsuba algorithm [11]. Despite the significant improvement gained by the Karatsuba algorithm, the complexity is still not optimal. Furthermore, the implementation of the Karatsuba algorithm is more burdensome due to its recursive nature. The known fastest multiplication algorithm, introduced by Sch¨ onhage and Strassen [21], performs multiplication in the frequency domain using the Fast Fourier Transform (FFT) [9] with complexity O(m log m log log m) for multiplication of m-bit integers or m-coefficient polynomials [10]. Unfortunately, the FFT based algorithm becomes efficient and useful in practice only for very large operands due to the overhead associated with the forward and inverse Fourier transform 1
operations. When the transformation computations are not considered, frequency domain polynomial multiplication (without modular reduction) has surprisingly low linear, i.e. O(m), complexity. Sadly, no efficient method for performing modular reduction is known to exist in the frequency domain, and therefore one needs to convert the product of a polynomial multiplication operation back to the time domain to perform modular reduction bearing significant overhead. Similar approaches on doing arithmetic over domains where multiplication has only linear complexity but reduction is difficult have been proposed with applications in cryptography [2],[3]. In [2], the authors propose a full implementation of RSA[20] in the Residue Number System (RNS) which does not require any conversions either from radix to RNS beforehand or from RNS to radix afterwards. In their implementation they use an RNS version of the Montgomery multiplication algorithm. In [3] an alternate polynomial representation is utilized where a polynomial is represented by its values at sufficiently many points and a Montgomery like modular multiplication algorithm is proposed for elliptic curve cryptography over the extension fields GF (pm ) for p > 2m. In this paper we propose a new algorithm for achieving modular reduction in the frequency domain. With our approach, the entire finite field multiplication, including modular reduction, can be carried out in the frequency domain. In many finite field applications a chain of arithmetic operations are performed rather than a solitary one. Using our method after an initial conversion step all intermediary operations may be computed in the frequency domain. Therefore, there will be no need for conversion before and after every single finite field multiplication except before the very first and after the very last ones. In the next section, we briefly give some background information on finite fields, the discrete Fourier transform (DFT) in finite fields and application of the DFT to multiplication. In Section 3, we propose a new algorithm, named DFT modular multiplication, for finite field modular multiplication of polynomials in the frequency domain. In Section 4, we present ideas for efficient implementation of DFT modular multiplication. Furthermore, in Section 5, we give a complexity analysis of the DFT modular multiplication algorithm, compare it with other efficient methods such as the Karatsuba algorithm, and show its relevance for elliptic curve cryptography. Finally, in Section 6 we present our conclusions.
2 2.1
Background Finite Fields and Polynomial Representation
A field with finite number of elements is called a finite field or Galois field, denoted by Fq or GF (q), where q stands for the number of elements in the field [14]. The number of elements in a finite field is always a prime or a prime power, i.e., q = p or q = pm ,where the prime number p is called the characteristic of the finite field. When q is a prime, i.e. q = p, the finite field GF (p) is called a prime field. The prime field GF (p) is the field of residue classes modulo p and its elements are represented by the integers in {0, 1, 2, . . . , p − 1}. When q is a prime power, i.e. q = pm , the finite field GF (pm ) is called an extension field. The extension field GF (pm ) is generated by using an mth degree irreducible polynomial over GF (p) and it is the field of residue classes modulo the irreducible field generating polynomial. Hence, in polynomial representation the elements of GF (pm ) are represented by polynomials of degree m − 1 with coefficients in GF (p).
2
2.2
Multiplication in Finite Fields Using the Schoolbook Method
Finite field polynomial multiplication is performed modulo Pm−1 thei field characteristic Pm−1and ithe field generating polynomial. For instance, let A(x) = a x and B(x) = i i=0 i=0 bi x , where ai , bi ∈ GF (q), be two polynomials representing two elements of the finite field GF (q m ) constructed with an mth degree field generating polynomial f (x) which is irreducible over GF (q). The most straightforward way of performing multiplication of two polynomials would be by the the schoolbook method of multiplication. One can perform the multiplication of A(x) and B(x) by the schoolbook method and find the product C(x) = A(x) · B(x) in two steps as follows: 1. Polynomial multiplication: 0
C (x) = A(x) · B(x) =
m−1 X m−1 X
ai bj xi+j
i=0 j=0
2. Modular reduction: C(x) = C 0 (x) (mod f (x))
Note above that C(x) is an (m − 1)st degree polynomial whose coefficients are elements of GF (q). As much as it is straightforward, schoolbook method of multiplication is inefficient and usually not the preferred method in finite field applications. In section 3, we will introduce an efficient method for finite field polynomial multiplication by using the Number Theoretic Transform [18], and compare it with other efficient algorithms such as the Karatsuba Multiplication [11] and a more recent method by Bajard et al. that utilizes an alternate polynomial representation and Lagrange’s theorem [3].
2.3
Number Theoretic Transform (NTT)
Number theoretic transform over a ring, also known as the discrete Fourier transform over a finite field, was introduced by Pollard [18]. For a finite field Fq and a generic sequence (a) of length d whose entries are from Fq , the forward DFT of (a) over Fq , denoted by (A), can be computed as Aj =
d−1 X
ai rij , 0 ≤ j ≤ d − 1 .
(1)
i=0
Here we refer to the elements of (a) and (A) by ai and Ai , respectively, for 0 ≤ i ≤ d − 1. Likewise, the inverse DFT of (A) over Fq can be computed as ai =
d−1 1 X · Aj r−ij , 0 ≤ i ≤ d − 1 . d
(2)
j=0
The above DFT computations over Fq are defined by utilizing a d-th primitive root of unity, denoted by r, from Fq or a finite extension of Fq . Definition 1 r is a primitive dth root of unity modulo n if rd = 1 (mod n) 3
and rd/t − 1 6= 0 (mod n) for any prime divisor t of d. Note that unlike the complex number r = ej2π/n used for the DFT computations in the complex number field, a finite field element r is used in finite fields. Also, the DFT over a finite field has no physical meaning as in the case of the DFT over the complex number field and serves only for the purpose of computing convolution efficiently as we will explain in section 2.3.2. However, in this paper we will still use the terms time domain and frequency domain for identifying a sequence and its number theoretic transform. Hence, for instance, we will refer to (a) and (A) as the time and frequency domain representations, respectively, of the same sequence. We would like to caution the reader that in a number theoretic transform the modulus q and the transform length d can not be chosen independently of each other. If the modulus is a composite number of the form q = pe11 pe22 pe33 . . . pekk , then the length of the transform, d, must divide gcd(p1 − 1, p2 − 1, p2 − 1, . . . , pk − 1). In the case of NTT over a finite field, i.e. when q is a prime p or a prime power pm , d must divide p − 1. 2.3.1
Mersenne Transform
Definition 2 A number of the form Mn = 2n −1, where n is a positive integer, is called a Mersenne number. If a Mersenne number is also a prime number, it is called a Mersenne prime. Mersenne primes are favourable choices as finite field characteristics due to their computational advantage in the modular reduction operation. Modular reduction by Mersenne primes can be achieved with simple addition and shift operations. A Number Theoretic Transform of special interest is Mersenne Transform, which is an NTT with arithmetic modulo a Mersenne number [19]. Mersenne Transform allows for very efficient forward and inverse transform operations for r = ±2. Multiplication of an n-bit number with integer powers of 2 modulo a Mersenne prime Mn can be performed by simple rotation of the bits of the n-bit number given as in the following properties. Property 1 For a positive integer i and an n-bit integer z ∈ GF (q), where q = 2n − 1 is a Mersenne prime, 2i · z mod q is equal to the n-bit number z rotated by i bits to the left. Property 2 For a positive integer i and an n-bit integer z ∈ GF (q), where q = 2n − 1 is a Mersenne prime, 2−i · z mod q is equal to the n-bit number z rotated by i bits to the right. Similarly, multiplication of an n-bit number with integer powers of −2 modulo a Mersenne prime Mn can be performed by simple rotation of the bits of the n-bit number, in addition to a negation if the power of −2 is odd. Hence, in a Mersenne transform with r = ±2 all of the multipllications by powers of r in the forward and inverse DFT computations can be achieved with simple bitwise rotations. In this case, for a transform length of d, the forward DFT computation can be achieved by only (d − 1)2 simple rotations and d(d − 1) additions/subtractions avoiding any multiplications. Likewise, the inverse DFT computation can be achieved by only (d − 1)2 simple rotations and d(d − 1) additions/subtractions, in addition to d constant multiplications by 1/d. For a more detailed complexity analysis of Mersenne transform, we refer the interested reader to [19]. Note that in a Mersenne transform the values of the sequence length d and the dth primitive root of unity r are dependent on each other. For the Mersenne transform with modulus q = Mn = 2n −1 and for r = ±2, the following equalities hold determining the relationship between d and r: 4
d=
n,
r=2.
2n ,
r = −2 .
For more information on several properties of Mersenne transform, the reader is referred to [19]. 2.3.2
Convolution Theorem and Multiplication in GF (q m ) Using NTT
A significant application of the Fourier transform is convolution. According to the convolution theorem, computing the convolution of two sequences in the time domain is equivalent to computing the componentwise parallel multiplication of their frequency domain representations. Convolution of any two d-element sequences (a) and (b) in the time domain results in another d-element sequence (c) and can be computed as follows: ci =
d−1 X
aj bi−j mod d , 0 ≤ i ≤ d − 1 .
(3)
j=0
The above convolution operation in the time domain is equivalent to the following computation in the frequency domain: Ci = Ai · Bi , 0 ≤ i ≤ d − 1 , (4) where (A), (B) and (C) denote the discrete Fourier transforms of (a), (b) and (c), respectively. Hence, convolution of two d-element sequences in the time domain, with complexity O(d2 ), is equivalent to simple pairwise multiplication of the DFTs of these sequences and has a surprisingly low O(d) complexity. Note in (3) that, the convolution is computed by making the sequences (a) and (b) periodic with d. This is equivalent to the cyclic convolution of the two sequences. Multiplication of two polynomials is basically the same as the acyclic (linear) convolution of the polynomial coefficients. We have seen that cyclic convolution can be performed very efficiently in the Fourier domain by pairwise coefficient multiplications. Hence, it will be wise to represent an element of GF (q m ), in polynomial representation an (m − 1)st degree polynomial with coefficients in GF (q), with at least a d = (2m − 1) element sequence by appending zeros at the end, so that the cyclic convolution of two such sequences will be equivalent to their acyclic convolution and give us their polynomial multiplication. We can form sequences by taking the ordered coefficients of polynomials. For instance, a(x) = a0 + a1 x + a2 x2 + . . . + am−1 xm−1 , an element of GF (q m ) in polynomial representation, can be interpreted as the following sequence after appending d − m zeros to the right: (a) = (a0 , a1 , a2 , . . . , am−1 , 0, 0, . . . , 0) .
(5)
For a(x), b(x) ∈ GF (q m ), and for d ≥ 2m − 1, the cyclic convolution of (a) and (b) yields a sequence (c) whose entries can be interpreted as the coefficients of a polynomial c(x) such that c(x) = a(x) · b(x). The computation of this cyclic convolution can be performed by simple pairwise coefficient multiplications in the discrete Fourier domain. The following straightforward algorithm realizes the polynomial multiplication c(x) = a(x) · b(x), without modular reduction by the field 5
generating polynomial, using the direct application of DFT. Algorithm for Polynomial Multiplication by Direct Application of DFT Input: a(x), b(x) ∈ GF (q m ) Output: c(x) = a(x) · b(x) Step 1. Interpret the coefficients of a(x) and b(x) as elements of the sequences (a) and (b), respectively, and append zeros to the right to make their lengths exactly d. Step 2. Convert (a) and (b) into their respective frequency domain representations (A) and (B) using the DFT operation of (1). Step 3. Multiply (A) and (B) together to compute (C) as in (4). Step 4. Convert (C) back to the time domain representation (c) using the inverse DFT operation of (2). Step 5. Interpret the first 2m − 1 coefficients of (c) as the coefficients of the product polynomial c(x). 2 The above method can be applied to multiplication of two polynomials of arbitrary degrees by only changing the condition d ≥ 2m − 1 to d > deg(a(x)) + deg(b(x)), where deg(p(x)) denotes the degree of p(x). Using this approach, polynomial multiplication is achieved with d(d − 1) additions and (d − 1)2 multiplications for the forward transform operation, d(d − 1) additions, (d − 1)2 multiplications and d constant multiplications for the inverse transform operation, and d multiplications for the multiplication in the frequency domain, totalling 2(d−1)2 +d multiplications, d constant multiplications and 2d(d − 1) additions in terms of coefficient arithmetic. However, as we mentioned in section 2.3.1, the same operations can be achieved with only 2(d − 1)2 bitwise rotations, d multiplications, d constant multiplications and 2d(d − 1) additions of coefficients when we use Mersenne transform. Note that, with this method the polynomial product c(x) = a(x) · b(x) is computed in the frequency domain but the final reduction by the field generating polynomial is not performed. One needs to convert (C) back to the time domain so that it can be reduced modulo the field generating polynomial and further multiplications can be performed on it using the same method. In the next section, we introduce a new method that performs both polynomial multiplication and modular reduction in the frequency domain and hence avoids costly back and forth conversions between the time and frequency domains.
3
Modular Multiplication in the Frequency Domain
In many finite field applications, a chain of arithmetic operations needs to be performed, rather than a solitary one. For example, in elliptic curve cryptography a scalar point product is computed by applying a chain of finite field additions, subtractions, multiplications, squarings and inversions on the input point coordinates [8]. The Montgomery residue representation has proven to be useful in this computation [16], [13], [12]. In using this method, first the operands are converted to their respective Montgomery residue representations, then utilizing Montgomery arithmetic the desired computation is implemented, and finally the result is converted back to the normal integer or polynomial representation. If there is a large number of operations performed in the Montgomery domain, due to the efficiency of the intervening computations, the forward and backward conversion operations become affordable. We introduce the same notion for the frequency domain. We want
6
to find an arithmetic operation in the frequency domain that is equivalent to ordinary modular or Montgomery multiplication in the time domain. For elliptic curve cryptography, we assume projective coordinates are used and all inversions are also performed by means of multiplications in the frequency domain. Thus, all finite field arithmetic operations, including modular reductions, may be performed in the frequency domain, avoiding any cost due to conversions between the time and the frequency domains. In the remainder of this section, we introduce the DFT modular multiplication algorithm which allows for both polynomial multiplication and modular reduction operations in the frequency domain. In section 3.1 we briefly present the notation we use in representing finite field elements in the frequency domain in the DFT modular multiplication algorithm and we present the algorithm in section 3.2.
3.1
Representation of Finite Field Elements in the DFT Modular Multiplication Algorithm
Let GF (q m ) denote a finite field extension generated by an irreducible polynomial f (x). For a(x), b(x) ∈ GF (q m ), both (m − 1)st degree polynomials with coefficients in GF (q) in polynomial representation, we let (a) and (b), respectively, denote their corresponding d-element time domain sequences. Note that the sequences (a) and (b) are extended to have length d by setting their higher d − m elements to 0. In this case, a(x) and b(x) have a maximum number of m coefficients each and, as explained in Sections 2 and 3, for the DFT operations we set d ≥ 2m − 1 and use r which is a d-th primitive root of unity in F = GF (q) or a finite extension of F. We denote the discrete Fourier transforms of (a) and (b) by (A) and (B), respectively, and use Ai and Bi , for 0 ≤ i ≤ d − 1, to denote their elements. Likewise, we represent the polynomial product c(x) = a(x) · b(x) with the d-element sequence (c). We let (C) represent the discrete Fourier transform of (c) and denote its elements by Ci , for 0 ≤ i ≤ d − 1. Finally, we let (f ) and (x) represent the d-element sequences for the irreducible field generating polynomial f (x) and the constant polynomial x ∈ GF (q m ), respectively, in the time domain. We denote the discrete Fourier transforms of (f ) and (x) by (F ) and (X), respectively, and let Fi and Xi , for 0 ≤ i ≤ d − 1, denote their elements. Hence, for instance, the constant polynomial x is represented by a d-element sequence as (x) = (0, 1, 0, 0, · · · , 0) , and the DFT of (x) is represented by the d-element sequence (X) given as (X) = (1, r, r2 , r3 , r4 , r5 , . . . , rd−1 ) . Note that we denote the elements of a sequence with the name of the sequence and a subscript for showing the location of the particular element in the sequence, e.g., X0 = 1 and Xd−1 = rd−1 .
3.2
The DFT Modular Multiplication Algorithm
In Table 1 we present the DFT modular multiplication algorithm which performs modular multiplication of a(x) and b(x), both elements of GF (q m ), in the discrete Fourier domain. We also give the complexity of each step of the algorithm in terms of the number of GF (q) multiplications and additions. The algorithm consists of two parts; multiplication (steps 1 and 2) and Montgomery reduction (steps 3 through 9). Multiplication is performed simply by pairwise multiplication of the elements of the frequency domain sequence representations of the two operands, i.e., Ai and Bi for 0 ≤ i ≤ d − 1. Multiplication of the two polynomials a(x) and b(x), each with degree at most (m − 1), results in another polynomial c(x) of degree at most 2m − 2. For performing further 7
multiplications over c(x) in the frequency domain, one needs to first reduce it modulo f (x) so that its time domain representation is of degree less than or equal to m − 1. Reduction is more complex and performed by Montgomery reduction in the frequency domain. The following algorithm achieves Montgomery reduction in our DFT modular multiplication algorithm. In the reduction process we use the normalized field generating polynomial f 0 (x) which equals f (x)/f0 mod q. Here, we let (f 0 ) denote the d-element time domain sequence representation for f 0 (x), with the sequence elements fi0 = fi · (f0−1 mod q) for 0 ≤ i ≤ d − 1. Hence, f 0 (x) is equivalent to f (x) but normalized to have its constant coefficient f 0 (0), or equivalently f00 in its time domain sequence representation (f 0 ), set to 1. Likewise, (F 0 ) denotes the DFT of (f 0 ) and its sequence elements are Fi0 = Fi · (f0−1 mod q), for 0 ≤ i ≤ d − 1. Note that, the reduction step of our algorithm results in the product c(x) of degree at most m − 1 by reducing the initial product modulo f 0 (x) which is equivalent to reducing it modulo f (x). However, it also adds an x−(m−1) factor to the otherwise correct result. DFT Montgomery Reduction Input: c(x) = a(x) · b(x) Output: c(x) = a(x) · b(x) · x−(m−1) mod f (x) Step 1. Compute −c0 , the negative of the first coefficient in the time domain sequence representation (c) of c(x), and put it in S (steps 4, 5, 6 and 7 in Table 1). Step 2. Make c0 zero by adding S · (f 0 ) to (c). Then shift the elements of new (c), which now has c0 = 0 and is a multiple of (x), by one position to the left via division by (x) (steps 8 and 9 in Table 1). Step 3. Repeat 1 and 2 for m − 1 times. Proof of Correctness: DFT Montgomery reduction is a direct adaptation of Montgomery reduction. We compute S such that (c(x) + S · f 0 (x)) is a multiple of x. We then divide (c(x) + S · f 0 (x)) by x and obtain a result which is congruent to c(x) · x−1 modulo f (x). By repeating this m − 1 times we obtain the final result which is congruent to c(x) · x−(m−1) modulo f (x). 2 The inputs of the DFT modular multiplication algorithm in Table 1 are the DFTs (A) and (B) of the d-element sequences (a) and (b) which represent the polynomials a(x) and b(x), respectively. The output of the algorithm is (C), the DFT of the sequence (c), where (c) represents the delement sequence for the polynomial c(x) = a(x)·b(x)·x−(m−1) mod f (x). The extra x−(m−1) factor shows that the DFT modular multiplication algorithm actually computes Montgomery products of polynomials [12] in the frequency domain. Hence, the input polynomials a(x) and b(x) may be viewed as the Montgomery residue representations of two polynomials u(x) and v(x) such that a(x) = u(x) · xm−1 mod f (x) and b(x) = v(x) · xm−1 mod f (x) . By the DFT modular multiplication the residue representation is kept intact, i.e.,
8
Step 1 Step 2 Step Step Step Step Step Step Step
3 4 5 6 7 8 9
MULTIPLICATION for i = 0 to d − 1 do Ci = Ai · Bi REDUCTION for j = 0 to m − 2 do S=0 for i = 0 to d − 1 do S = S + Ci S = −S/d for i = 0 to d − 1 do Ci = (Ci + Fi0 · S) · Xi−1
#Multiplication – d
#Add./Subtract. – –
– – – – m−1 – 2 · (m − 1) · d
– – – (m − 1) · d – – (m − 1) · d
Table 1: DFT modular multiplication algorithm for modular multiplication in GF (q m ) in the discrete Fourier domain
a(x) · b(x) · x−(m−1) mod f (x) = (u(x) · v(x)) · xm−1 mod f (x) which allows for further computation in the frequency domain. At first glance the complexity of the DFT modular multiplication algorithm appears to be even worse than that of the schoolbook multiplication method. For d ≈ 2m, modular multiplication in GF (q m ) with the DFT modular multiplication method requires 4m2 − m − 1 multiplications and 4m2 − 4m additions in the ground field GF (q), while the schoolbook method requires only m2 multiplications and (m − 1)2 additions ignoring the cost of modular reduction. Nevertheless, this complexity may be improved dramatically by using special values for q, r, d and the irreducible field generating polynomial f (x) as we will explain in the next section.
4
Efficient Implementation of DFT Modular Multiplication over Finite Fields Constructed with Special Parameters
In this section, we investigate the existence of finite fields of the form GF (q m ), where q is a Mersenne prime, and efficient parameters such as irreducible field generating polynomials of special form, sequence lengths satisfying d ≥ 2m − 1 and d ≈ 2m, and d-th primitive roots of unity r of special form exist for efficient implementation of DFT modular multiplication. We would like to caution the reader that all these parameters are dependent on each other and can not be chosen independently.
4.1
Sequence Length d
In section 2.3.2 we have seen that a necessary condition for performing polynomial multiplication in the frequency domain is that d ≥ 2m − 1. Using the smallest possible sequence length d will lead to smallest number of arithmetic operations in the computation of DFT modular multiplication (see Table 1). Optimally, d = 2m − 1 will lead to the least number of arithmetic operations. Furthermore, note that multiplications by −d−1 in step 7 of the DFT modular multiplication algorithm are constant multiplications, and hence can be computed by a sequence of additions and shift operations. For an n-bit field characteristic q, these constant multiplications can be achieved with n/2 shifts and n/2 additions, on average. And naturally, for lower Hamming weight 9
−d−1 values, the constant multiplications by −d−1 can be performed more efficiently with lower number of shifts and additions. This would be very helpful especially for implementations on processors without a multiplication instruction where multiplications are performed with shift and add instructions. However, again we would like to note that the parameter d can not be selected independently and its value depends on other parameters such as the field characteristic q and the dth primitive root of unity r.
4.2
Mersenne Primes as Finite Field Characteristics
As mentioned in section 2.3.1, the use of a Mersenne number as the modulus and selection of r = ±2 allow for very efficient DFT computations in Mersenne transform. We will utilize this fact in our DFT modular multiplication algorithm. In a finite field with a Mersenne prime characteristic, a multiplication by ri , for r = ±2 and i an integer, can be achieved by a simple bitwise rotation (in addition to a negation when r = −2 and i is odd) which is inexpensive (see section 2.3.1). In our DFT modular multiplication algorithm of Table 1, this property may be exploited if we choose q to be a Mersenne prime. In step 9 of the algorithm, Xi−1 = r−i and hence multiplications by Xi−1 can be computed with simple bitwise rotations. Multiplications by Fi0 can also be avoided in a similar fashion. For instance, for f (x) = xm + rs0 with s0 an integer, Fi0 = rmi−s0 mod d + 1 ,
(6)
and hence for r = ±2 multiplications by Fi0 can be performed with only one bitwise rotation and one addition/subtraction.
4.3
Binomials or Trinomials as Field Generating Polynomials
The efficiency of step 9 of the DFT modular multiplication algorithm can be improved significantly by using efficient parameters for r and the field generating polynomial f (x). The critical operations in step 9 are the two multiplications, which are multiplications of Fi0 by S and Ci + Fi0 · S by Xi−1 . The multipliers Fi0 are constants and their values depend on the irreducible polynomial f (x), as well as on the value of r. Similarly, the multipliers Xi−1 are constants whose values depend on the value of r. When r = 2, Xi−1 becomes r−i = 2−i = 2d−i and hence when the finite field characteristic is a Mersenne prime multiplications by Xi−1 , as in step 9 of the DFT modular multiplication algorithm, become simple (d − i)-bit left-rotations (or equivalently i-bit right-rotations), which have negligible cost compared to regular multiplication. Similarly, when r = −2, multiplications by Xi−1 = r−i = (−2)−i = (−2)d−i become simple (d − i)-bit left-rotations (or equivalently i-bit right-rotations) and a negation depending on whether the power of −2 is odd. The multiplications by Fi0 in step 9 of the DFT modular multiplication algorithm can be simplified significantly by using a low Hamming weight field generating polynomial for f (x). We can guarantee that the elements of the sequence (F ) are in the form of summations of powers of r by selecting the coefficients of f (x) to be 1 or, for that matter, any integer power of r. When f0 and all other coefficients of f (x) are powers of r, so are the coefficients of f 0 (x) = f (x)/f0 . In this case, the elements of the sequence (F 0 ) will also be in the form of summations of powers of r. Hence, for r = ±2 and in a finite field with a Mersenne prime characteristic, multiplications by Fi0 can be performed with simple bitwise rotations and additions/subtractions. For clarification, we make the following illustration. For an irreducible field generating polynomial of the form f (x) = xm + rsm−1 xm−1 + rsm−2 xm−2 + rsm−3 xm−3 + · · · + rs2 x2 + rs1 x + rs0 , 10
whose coefficients are all integer powers of r, the d-element time domain sequence representation is (f ) = (rs0 , rs1 , rs2 , rs3 , · · · , rsm−2 , rsm−1 , 1, 0, 0, · · · , 0) and the DFT of (f ) is (F ) = (F0 , F1 , F2 , F3 , · · · , Fd−1 ) , where Fi = rmi + r(m−1)i+sm−1 + r(m−2)i+sm−2 + r(m−3)i+sm−3 + · · · + r2i+s2 + ri+s1 + rs0 . Likewise, for the normalized field generating polynomial f 0 (x) = f (x)/f0 = r−s0 xm + rsm−1 −s0 xm−1 + rsm−2 −s0 xm−2 + rsm−3 −s0 xm−3 + · · · + rs1 −s0 x + 1 , the d-element time domain sequence representation is (f 0 ) = (1, rs1 −s0 , rs2 −s0 , rs3 −s0 , · · · , rsm−2 −s0 , rsm−1 −s0 , r−s0 , 0, 0, · · · , 0) and the DFT of (f 0 ) is 0 (F 0 ) = (F00 , F10 , F20 , F30 , · · · , Fd−1 ),
where Fi0 = Fi /f0 = rmi−s0 + r(m−1)i+sm−1 −s0 + r(m−2)i+sm−2 −s0 + r(m−3)i+sm−3 −s0 + · · · + ri+s1 −s0 + 1 . Note that all elements of the sequence (F 0 ), namely Fi0 for 0 ≤ i ≤ d − 1, are in the form of summations of powers of r. Therefore, when a finite field with a Mersenne prime characteristic is used and for r = ±2, multiplications by Fi0 can be performed with simple bitwise rotations and additions/subtractions. Clearly, selecting low Hamming weight polynomials for f (x) may reduce the complexity of multiplications by Fi0 by means of reducing the number of required additions/subtractions, as we shall see better next. In Table 2 we list the DFTs for some simple polynomials. Our intention is to determine the type of irreducible polynomial that has a simpler DFT and therefore works best for our purpose. Since the DFT is a linear operation, the DFT of the sum of two polynomials equals the sum of their individual DFTs. Similarly, the DFT of a constant multiple of a polynomial equals the same constant multiple of the DFT of that polynomial. Hence, the DFT of any polynomial is a summation of constant multiples of sequences such as those given in Table 2. f (x) r x xm
(F ) ( r ( 1 ( 1
, , ,
r r rm
, , ,
r r2 r2m
, , ,
r r3 r3m
, , ,
r r4 r4m
, , ,
... ... ...
, , ,
r ) rd−1 ) r(d−1)m )
Table 2: Frequency domain sequence representations for some simple polynomials Naturally, an irreducible binomial of the form f (x) = xm ± rs0 with an integer s0 , or trinomial 0 of the form f (x) = xm ± rsm0 xm ± rs0 with m0 < m and integers sm0 and s0 , will have a simpler frequency domain sequence representation and multiplications by F 0 , as in step 9 of the DFT 11
modular multiplication algorithm in Table 1, will be easier to perform. For instance, provided that f (x) = xm + rs0 we will have Fi = rmi + rs0 and Fi0 = f0−1 Fi = rmi−s0 + 1 . And for r = ±2 and a Mersenne prime q, multiplication of a number in GF (q) by Fi0 , of the form rmi−s0 + 1, can easily be performed with only a bitwise rotation and an addition/subtraction. 0 Likewise, for the trinomial f (x) = xm + rsm0 xm + rs0 , 0
Fi0 = rmi−s0 + rm i+sm0 −s0 + 1 , and multiplications by Fi0 can be performed with only two bitwise rotations and two additions/subtractions. Also, note that for any integer s the equality rs = rs mod d holds, since r is a dth primitive root of unity in GF (q), a condition satisfied for the application of the DFT, and therefore all required bitwise rotations due to multiplications by powers of r will be for less than d bits. Moreover, since r is a dth primitive root of unity, all numbers of the form rs , for any integer s, will have only d distinct values. Hence, in step 9 of the DFT modular multiplication algorithm, all but d of the multiplications of S by rk , due to multiplications by Fi0 and Xi−1 , may be avoided by precomputing the d distinct values of the rk · S products, for the value of S computed at step 7, and using these precomputed values as needed in the computations of Ci = (Ci + Fi0 · S) · Xi−1 .
4.4
Existence of Efficient Parameters
In Table 3 we give a list of parameters such as Mersenne prime Mn , prime extension degree m, sequence length d, r = ±2 and equivalent binary field size for which efficient modular multiplication in GF (q m ) using the DFT modular multiplication algorithm is possible for elliptic curve cryptography. Note that between d and n we always have the relationships d = n or d = 2n for r = 2 or r = −2, respectively. Also note that in our list we include only those cases where the field extension degree m is prime, which is a security requirement for elliptic curve cryptography over extension degrees. For each case listed in Table 3, in Appendix we also provide a list of irreducible binomials or trinomials as field generating polynomials that would allow for efficient implementation of DFT modular multiplication. For some of the cases we were able to find efficient irreducible binomials and included them, while for other cases we were not able to find such binomials and included trinomials instead. However, we would like to note that for the computationally more desirable cases of m = d/2 we were always able to find efficient irreducible binomials.
5
Complexity of the DFT Modular Multiplication Algorithm
In this section, we present the complexity of the DFT modular multiplication algorithm for a practical set of parameters and compare it with three other efficient multiplication approaches, namely the direct application of DFT multiplication with time domain modular reduction (see section 2.3.2), algorithm by Bajard et al. [3] and Karatsuba algorithm [11]. In our complexity 12
n 13 17 17 17 19 19 19 19 31 31
q = Mn 8191 131071 131071 131071 524287 524287 524287 524287 2147483647 2147483647
m 13 11 13 17 11 13 17 19 11 13
d 26 34 34 34 38 38 38 38 31 31
r −2 −2 −2 −2 −2 −2 −2 −2 2 2
equivalent binary field size ∼ 2169 ∼ 2187 ∼ 2221 ∼ 2289 ∼ 2209 ∼ 2247 ∼ 2323 ∼ 2361 ∼ 2341 ∼ 2403
Table 3: List of Mn , m, d and r = ±2 values that would make efficient DFT modular multiplication in GF (q m ) possible for elliptic curve cryptography over finite fields of size 150 to 500 bits.
1 2 3 4 5 6 7 8 9
MULTIPLICATION for i = 0 to d − 1 do Ci = Ai · Bi REDUCTION for j = 0 to m − 2 do S=0 for i = 0 to d − 1 do S = S + Ci S = −S/d for i = 0 to d − 1 do Ci = (Ci + Fi0 · S) · Xi−1
#Mult. – d ≈ 2m
#Const. Mult. – –
#Add./Subtract. – –
#Rotation – –
– – – –
–
– – – –
– – –
m−1 – –
– – – (m − 1) · d ≈ 2m(m − 1) – – 2(m − 1) · d ≈ 4m(m − 1)
– –
– – 2 · (m − 1) · d ≈ 4m(m − 1)
Table 4: DFT modular multiplication algorithm in GF (q m ) and its complexity for a Mersenne prime q, an irreducible field generating binomial of the form f (x) = xm + rs0 with an integer s0 , r = ±2 and d ≈ 2m.
analysis we assume the use of efficient parameters such as Mersenne prime finite field characteristics as q, irreducible field generating binomials of the form f (x) = xm + 2s0 with a positive integer s0 , d-th primitive root of unity r = ±2 and sequence length d ≈ 2m. The field parameters we use, such as low Hamming weight field generating polynomial and Mersenne prime field characteristics, would lead to efficient implementation of multiplication for all four methods. Therefore, for the selected parameters, we can safely assume that our comparisons are fair. In our complexity derivations for the DFT modular multiplication we utilize all efficient implementation ideas mentioned in Section 4 and obtain the complexity of multiplication in GF (q m ) in terms of the number of GF (q) operations as shown in Table 4. Note that, for performing DFT modular multiplication in GF (q m ) only a linear number of GF (q) multiplications (step 2) and constant multiplications (step 7) are needed, in addition to some additions/subtractions (steps 6, 9) and rotation operations (step 9) in GF (q) which are relatively inexpensive. In Table 5, we present the complexities of GF (q m ) multiplication in terms of GF (q) operations for all four approaches. The complexity presented for Karatsuba algorithm is only an approximate one for the case when m is a power of two. For the best complexities of Karatsuba algorithm for 13
#Multiplication #Const. Mult. #Add./Subtr. #Rotation
DFT Multiplication (direct approach) d ≈ 2m d ≈ 2m 2d(d − 1) ≈ 8m2 − 4m 2(d − 1)2 ≈ 8m2 − 8m + 2
Bajard et al.
Karatsuba
3m 2m + 2m
≈ mlog2 3 –
2m2 − m
≈ 6mlog2 3 − 7m + 1
–
m−1
2
DFT Modular Multiplication d ≈ 2m m−1 3(m − 1)d ≈ 6m2 − 6m 2(m − 1)d ≈ 4m2 − 4m
Table 5: Complexities of modular multiplication in GF (q m ) where the field generating polynomial is f (x) = xm + 2s0 for an integer s0 , q is a Mersenne prime and d ≈ 2m
different extension degrees, the reader is referred to [22]. Note that the number of ground field multiplications required by Karatsuba algorithm is subquadratic, whereas all other approaches require only linear number of ground field multiplications. Furthermore, although the algorithm by Bajard et al. has only linear complexity in terms of the number of multiplications, it still requires performing a quadratic number of constant multiplications in GF (q). Whereas the number of constant multiplications required by DFT modular multiplication is only linear. In the number of additions and subtractions, Karatsuba algorithm has subquadratic complexity whereas the complexities of the other three methods are quadratic. Finally, multiplication with the direct DFT approach and DFT modular multiplication algorithm require a quadratic number of simple bitwise rotations, whereas Karatsuba multiplication requires only a linear number of bitwise rotations and the method by Bajard et al. does not require any. Clearly, the complexity of the DFT modular multiplication algorithm is better than the direct DFT approach. Moreover, since DFT modular multiplication requires significantly less number of complex operations such as multiplication and constant multiplication, its overall performance appears to be better than the algorithm by Bajard et al. and the Karatsuba algorithm, especially for computational platforms where multiplication is expensive compared to other operations such as addition, subtraction or bitwise rotation. In many modern microprocessors speed of a multiplication instruction is comparable to or only a few times slower than an addition instruction. However this is not the case especially for some microprocessors where there are no special multiplication instructions and multiplications are performed with a sequence of shifts and additions. For instance, ATMEL 90LS8535 [1] is a microprocessor that is widely used in embedded applications, yet it has no multiplication instructions. Likewise, in hardware complexity of multiplication is significantly more than complexity of addition. For instance, a straightforward implementation of an n-bit multiplication can be achieved via n additions and n shift operations. Therefore, in serialized implementation the complexity of an n-bit multiplication may be assumed to be roughly equal to the complexity of n n-bit additions and n shift operations. Using the same approach, we may assume that multiplication by a constant n-bit number can be performed with n/2 shifts and n/2 additions on average, and hence has complexity equal to n/2 addition and n/2 shift operations. Under these assumptions, Table 5 gives the complexities of modular multiplication in GF (q 13 ) for all four methods where q = 213 − 1 and GF (q 13 ) is constructed using the irreducible binomial f (x) = x13 − 2. The table also includes the total number of clock cycles that a single multiplication with these methods takes, assuming addition/subtraction and rotation operations all take a single clock cycle. Note that this finite field has size ∼ 2169 and can be considered suitable for elliptic curve cryptography. We caution the reader that all values in Table 5 are obtained directly from
14
#Multiplication #Const. Mult. #Add./Subtr. #Rotation #Total Clock Cycles
DFT Multiplication (direct approach) 26 26 1300 1250 3564
Bajard et al.
Karatsuba
39 364 325 – 6071
91 – 390 12 2768
DFT Modular Multiplication 26 12 936 624 2392
Table 6: Complexities of modular multiplication in GF (q 13 ) where the field generating polynomial is f (x) = x13 − 2, q = 213 − 1, d = 26 and r = −2
Table 5, except the complexity of Karatsuba multiplication for GF (q 13 ) which we draw from [22] . Also, note that for all four methods we assume the average case for the complexities of constant multiplications regardless of the values of the particular parameters used. Hence, as demonstrated with the examples above, finite field polynomial multiplication may be achieved more efficiently with DFT modular multiplication than with other efficient methods such as the Karatsuba algorithm in computational platforms where multiplication operation is significantly more expensive than simple operations such as addition, subtraction and bitwise rotation.
6
Conclusion
We introduced the DFT modular multiplication algorithm which performs modular multiplication in the frequency domain using Montgomery reduction. The FFT based multiplication technique, which performs multiplication in the frequency domain, is known to be very efficient. Nevertheless, due to the lack of a frequency domain reduction algorithm, for performing modular reductions conversion to the time domain is needed, which adds a significant overhead and makes the FFT based multiplication algorithms impractical for small operands. In this work, by allowing for modular reductions in the frequency domain, our algorithm avoids back and forth conversions between the frequency and the time domains. We have shown that, with our method, in computational platforms where multiplication is expensive frequency domain modular multiplication for finite fields becomes practical and can be achieved more efficiently than multiplication in the time domain for finite field sizes relevant to elliptic curve cryptography.
References [1] Atmel AVR 8bit RISC processor. http://www.atmel.com/dyn/resources/prod documents/1041S.PDF. [2] J.-C. Bajard and L. Imbert. A Full RNS Implementation of RSA. IEEE Transactions on Computers, 53(6):769–774, June 2004. [3] J.-C. Bajard, L. Imbert, C. N`egre, and T. Plantard. Efficient Multiplication in GF (pk ) for Elliptic Curve Cryptography. In Proceedings of the 16th IEEE Symposium on Computer Arithmetic (ARITH’03). IEEE, 2003. [4] E. R. Berlekamp. Algebraic Coding Theory. McGraw-Hill, New York, New York, USA, 1968. 15
[5] R. E. Blahut. Theory and Practice of Error Control Codes. Addison-Wesley, Reading, Massachusetts, USA, 1983. [6] R. E. Blahut. Fast Algorithms for Digital Signal Processing. Addison-Wesley, Reading, Massachusetts, USA, 1985. [7] I.F. Blake, X.H. Gao, R.C. Mullin, S.A. Vanstone, and T. Yaghgoobin. Applications of Finite Fields. Kluwer Academic, 1999. [8] I.F. Blake, G. Seroussi, and N. Smart. Elliptic Curves in Cryptography. Cambridge University Press, London Mathematical Society Lecture Notes Series 265, 1999. [9] J. Cooley and J. Tukey. An Algorithm for the Machine Calculation of Complex Fourier Series. Mathematics of Computation, 19:297–301, 1965. [10] R. Crandall and C. Pomerance. Prime Numbers. Springer-Verlag, New York, NY, USA, 2001. [11] A. Karatsuba and Y. Ofman. Multiplication of Multidigit Numbers on Automata. Sov. Phys. Dokl. (English translation), 7(7):595–596, 1963. [12] C ¸ . K Ko¸c and T. Acar. Montgomery Multplication in GF (2k ). Design, Codes, and Cryptography, 14(1):57–69, 1998. ¸ . K Ko¸c, T. Acar, and B. Kaliski. Analyzing and Comparing Montgomery Multiplication [13] C Algorithms. IEEE Micro, pages 26–33, June 1996. [14] R. Lidl and H. Niederreiter. Finite Fields, volume 20 of Encyclopedia of Mathematics and its Applications. Addison-Wesley, Reading, Massachusetts, USA, 1983. [15] R. J. McEliece. Finite Fields for Computer Scientists and Engineers. Kluwer Academic Publishers, 2nd edition, 1989. [16] P. L. Montgomery. Modular Multiplication without Trial Division. Mathematics of Computation, 44(170):519–521, April 1985. [17] C. Paar. Efficient VLSI Architectures for Bit-Parallel Computation in Galois Fields. PhD thesis, (Engl. transl.), Institute for Experimental Mathematics, University of Essen, Essen, Germany, June 1994. ISBN 3–18–332810–0. [18] J. M. Pollard. The Fast Fourier Transform in a Finite Field. Mathematics of Computation, 25:365–374, 1971. [19] C. M. Rader. Discrete Convolutions via Mersenne Transforms. IEEE Transactions on Computers, C-21(12):1269–1273, December 1972. [20] R. L. Rivest, A. Shamir, and L. Adleman. A Method for Obtaining Digital Signatures and Public-Key Cryptosystems. Communications of the ACM, 21(2):120–126, February 1978. [21] A. Sch¨ onhage and V. Strassen. Schnelle Multiplikation großer Zahlen. Computing, 7:281–292, 1971. [22] A. Weimerskirch and C. Paar. Generalizations of the Karatsuba Algorithm for Efficient Implementations. Technical report, Department of Electrical Engineering and Information Sciences, Ruhr-Universit¨at Bochum, Germany. 16
Appendix x13 − 2 x13 − 22 x13 − 23 x13 − 24 x13 − 25 x13 − 26
x13 − 27 x13 − 28 x13 − 29 x13 − 210 x13 − 211 x13 − 212
x13 + 2 x13 + 22 x13 + 23 x13 + 24 x13 + 25 x13 + 26
x13 + 27 x13 + 28 x13 + 29 x13 + 210 x13 + 211 x13 + 212
Table 7: Irreducible binomials of the form x13 ±2s0 for 0 ≤ s0 ≤ 12 that can be used for construction of the finite field GF (q 13 ) where q = 213 − 1.
x11 + 28 x3 + 1 x11 + 211 x3 + 22 x11 + 22 x3 + 24 x11 + x3 + 26 x11 + 23 x3 + 28 x11 + 26 x3 + 210 x11 + 29 x3 + 212 x11 + x3 + 214 x11 + 23 x3 + 216
x11 + 213 x3 + 1 x11 + 216 x3 + 22 x11 + 214 x3 + 24 x11 + 25 x3 + 26 x11 + 28 x3 + 28 x11 + 211 x3 + 210 x11 + 214 x3 + 212 x11 + 212 x3 + 214 x11 + 215 x3 + 216
x11 + 2x3 + 2 x11 + 24 x3 + 23 x11 + 27 x3 + 25 x11 + 210 x3 + 27 x11 + 2x3 + 29 x11 + 24 x3 + 211 x11 + 22 x3 + 213 x11 + 25 x3 + 215
x11 + 26 x3 + 2 x11 + 29 x3 + 23 x11 + 212 x3 + 25 x11 + 215 x3 + 27 x11 + 213 x3 + 29 x11 + 216 x3 + 211 x11 + 27 x3 + 213 x11 + 210 x3 + 215
Table 8: Irreducible trinomials of the form x11 + 2s3 x3 + 2s0 for 0 ≤ s3 , s0 ≤ 16 that can be used for construction of the finite field GF (q 11 ) where q = 217 − 1.
17
x13 + 2x + 1 x13 + 212 x + 22 x13 + 26 x + 24 x13 + x + 26 x13 + 211 x + 28 x13 + 25 x + 210 x13 + x + 212 x13 + 210 x + 214 x13 + 24 x + 216
x13 + 22 x + 1 x13 + 213 x + 22 x13 + 27 x + 24 x13 + 2x + 26 x13 + 212 x + 28 x13 + 26 x + 210 x13 + 216 x + 212 x13 + 211 x + 214 x13 + 25 x + 216
x13 + 215 x + 2 x13 + 29 x + 23 x13 + 23 x + 25 x13 + 214 x + 27 x13 + 28 x + 29 x13 + 22 x + 211 x13 + 213 x + 213 x13 + 27 x + 215
x13 + 216 x + 2 x13 + 210 x + 23 x13 + 24 x + 25 x13 + 215 x + 27 x13 + 29 x + 29 x13 + 23 x + 211 x13 + 214 x + 213 x13 + 28 x + 215
Table 9: Irreducible trinomials of the form x13 + 2s1 x + 2s0 for 0 ≤ s1 , s0 ≤ 16 that can be used for construction of the finite field GF (q 13 ) where q = 217 − 1.
x17 − 2 x17 − 22 x17 − 23 x17 − 24 x17 − 25 x17 − 26 x17 − 27 x17 − 28
x17 − 29 x17 − 210 x17 − 211 x17 − 212 x17 − 213 x17 − 214 x17 − 215 x17 − 216
x17 + 2 x17 + 22 x17 + 23 x17 + 24 x17 + 25 x17 + 26 x17 + 27 x17 + 28
x17 + 29 x17 + 210 x17 + 211 x17 + 212 x17 + 213 x17 + 214 x17 + 215 x17 + 216
Table 10: Irreducible binomials of the form x17 ± 2s0 for 0 ≤ s0 ≤ 16 that can be used for construction of the finite field GF (q 17 ) where q = 217 − 1.
x11 + 29 x + 1 x11 + 2x + 22 x11 + 24 x + 24 x11 + 211 x + 26 x11 + 23 x + 28 x11 + 26 x + 210 x11 + 213 x + 212 x11 + 2x + 214 x11 + 28 x + 216 x11 + x + 218
x11 + 213 x + 1 x11 + 216 x + 22 x11 + 28 x + 24 x11 + 215 x + 26 x11 + 218 x + 28 x11 + 210 x + 210 x11 + 217 x + 212 x11 + 25 x + 214 x11 + 212 x + 216 x11 + 215 x + 218
x11 + 23 x + 2 x11 + 210 x + 23 x11 + 22 x + 25 x11 + 25 x + 27 x11 + 212 x + 29 x11 + x + 211 x11 + 27 x + 213 x11 + 214 x + 215 x11 + 22 x + 217
x11 + 27 x + 2 x11 + 214 x + 23 x11 + 217 x + 25 x11 + 29 x + 27 x11 + 216 x + 29 x11 + 24 x + 211 x11 + 211 x + 213 x11 + 218 x + 215 x11 + 26 x + 217
Table 11: Irreducible trinomials of the form x11 + 2s1 x + 2s0 for 0 ≤ s1 , s0 ≤ 18 that can be used for construction of the finite field GF (q 11 ) where q = 219 − 1.
18
x13 + 24 x + 1 x13 + 215 x + 24 x13 + 27 x + 28 x13 + 218 x + 212 x13 + 210 x + 216
x13 + 22 x + 2 x13 + 213 x + 25 x13 + 25 x + 29 x13 + 216 x + 213 x13 + 28 x + 217
x13 + x + 22 x13 + 211 x + 26 x13 + 23 x + 210 x13 + 214 x + 214 x13 + 26 x + 218
x13 + 217 x + 23 x13 + 29 x + 27 x13 + 2x + 211 x13 + 212 x + 215
Table 12: Irreducible trinomials of the form x13 + 2s1 x + 2s0 for 0 ≤ s1 , s0 ≤ 18 that can be used for construction of the finite field GF (q 13 ) where q = 219 − 1.
x17 + 25 x + 1 x17 + 211 x + 24 x17 + 217 x + 28 x17 + 24 x + 212 x17 + 210 x + 216
x17 + 216 x + 2 x17 + 23 x + 25 x17 + 29 x + 29 x17 + 215 x + 213 x17 + 22 x + 217
x17 + 28 x + 22 x17 + 214 x + 26 x17 + 2x + 210 x17 + 27 x + 214 x17 + 213 x + 218
x17 + x + 23 x17 + 26 x + 27 x17 + 212 x + 211 x17 + 218 x + 215
Table 13: Irreducible trinomials of the form x17 + 2s1 x + 2s0 for 0 ≤ s1 , s0 ≤ 18 that can be used for construction of the finite field GF (q 17 ) where q = 219 − 1.
x19 − 2 x19 − 22 x19 − 23 x19 − 24 x19 − 25 x19 − 26 x19 − 27 x19 − 28 x19 − 29
x19 − 210 x19 − 211 x19 − 212 x19 − 213 x19 − 214 x19 − 215 x19 − 216 x19 − 217 x19 − 218
x19 + 2 x19 + 22 x19 + 23 x19 + 24 x19 + 25 x19 + 26 x19 + 27 x19 + 28 x19 + 29
x19 + 210 x19 + 211 x19 + 212 x19 + 213 x19 + 214 x19 + 215 x19 + 216 x19 + 217 x19 + 218
Table 14: Irreducible binomials of the form x19 ± 2s0 for 0 ≤ s0 ≤ 18 that can be used for construction of the finite field GF (q 19 ) where q = 219 − 1.
19
x11 + 220 x + 1 x11 + 219 x + 22 x11 + 218 x + 24 x11 + 217 x + 26 x11 + 216 x + 28 x11 + 215 x + 210 x11 + 214 x + 212 x11 + 213 x + 214 x11 + 212 x + 216 x11 + 211 x + 218 x11 + 210 x + 220 x11 + 29 x + 222 x11 + 28 x + 224 x11 + 27 x + 226 x11 + 26 x + 228 x11 + 25 x + 230
x11 + 229 x + 1 x11 + 228 x + 22 x11 + 227 x + 24 x11 + 226 x + 26 x11 + 225 x + 28 x11 + 224 x + 210 x11 + 223 x + 212 x11 + 222 x + 214 x11 + 221 x + 216 x11 + 220 x + 218 x11 + 219 x + 220 x11 + 218 x + 222 x11 + 217 x + 224 x11 + 216 x + 226 x11 + 215 x + 228 x11 + 214 x + 230
x11 + 24 x + 2 x11 + 23 x + 23 x11 + 22 x + 25 x11 + 2x + 27 x11 + x + 29 x11 + 28 x + 211 x11 + 27 x + 213 x11 + 26 x + 215 x11 + 25 x + 217 x11 + 24 x + 219 x11 + 23 x + 221 x11 + 22 x + 223 x11 + 2x + 225 x11 + x + 227 x11 + 221 x + 229
x11 + 213 x + 2 x11 + 212 x + 23 x11 + 211 x + 25 x11 + 210 x + 27 x11 + 29 x + 29 x11 + 230 x + 211 x11 + 229 x + 213 x11 + 228 x + 215 x11 + 227 x + 217 x11 + 226 x + 219 x11 + 225 x + 221 x11 + 224 x + 223 x11 + 223 x + 225 x11 + 222 x + 227 x11 + 230 x + 229
Table 15: Irreducible trinomials of the form x11 + 2s1 x + 2s0 for 0 ≤ s1 , s0 ≤ 30 that can be used for construction of the finite field GF (q 11 ) where q = 231 − 1.
20
x13 + 28 x + 1 x13 + 2x + 2 x13 + 22 x + 22 x13 + 26 x + 23 x13 + 22 x + 24 x13 + x + 25 x13 + 24 x + 26 x13 + x + 27 x13 + 213 x + 28 x13 + 22 x + 29 x13 + 27 x + 210 x13 + 211 x + 211 x13 + x + 212 x13 + 25 x + 213 x13 + 29 x + 214 x13 + 22 x + 215 x13 + 23 x + 216 x13 + 27 x + 217 x13 + x + 218 x13 + 2x + 219 x13 + 25 x + 220 x13 + 2x + 221 x13 + 214 x + 222 x13 + 23 x + 223 x13 + 28 x + 224 x13 + 212 x + 225 x13 + 2x + 226 x13 + 26 x + 227 x13 + 210 x + 228 x13 + 23 x + 229 x13 + 24 x + 230
x13 + 212 x + 1 x13 + 24 x + 2 x13 + 217 x + 22 x13 + 210 x + 23 x13 + 211 x + 24 x13 + 215 x + 25 x13 + 28 x + 26 x13 + 29 x + 27 x13 + 217 x + 28 x13 + 26 x + 29 x13 + 222 x + 210 x13 + 215 x + 211 x13 + 24 x + 212 x13 + 220 x + 213 x13 + 213 x + 214 x13 + 25 x + 215 x13 + 218 x + 216 x13 + 211 x + 217 x13 + 23 x + 218 x13 + 216 x + 219 x13 + 29 x + 220 x13 + 210 x + 221 x13 + 218 x + 222 x13 + 27 x + 223 x13 + 223 x + 224 x13 + 216 x + 225 x13 + 25 x + 226 x13 + 221 x + 227 x13 + 214 x + 228 x13 + 26 x + 229 x13 + 219 x + 230
x13 + 215 x + 1 x13 + 213 x + 2 x13 + 221 x + 22 x13 + 213 x + 23 x13 + 226 x + 24 x13 + 219 x + 25 x13 + 211 x + 26 x13 + 224 x + 27 x13 + 220 x + 28 x13 + 29 x + 29 x13 + 226 x + 210 x13 + 218 x + 211 x13 + 27 x + 212 x13 + 224 x + 213 x13 + 216 x + 214 x13 + 214 x + 215 x13 + 222 x + 216 x13 + 214 x + 217 x13 + 212 x + 218 x13 + 220 x + 219 x13 + 212 x + 220 x13 + 225 x + 221 x13 + 221 x + 222 x13 + 210 x + 223 x13 + 227 x + 224 x13 + 219 x + 225 x13 + 28 x + 226 x13 + 225 x + 227 x13 + 217 x + 228 x13 + 215 x + 229 x13 + 223 x + 230
x13 + 224 x + 1 x13 + 228 x + 2 x13 + 224 x + 22 x13 + 222 x + 23 x13 + 230 x + 24 x13 + 222 x + 25 x13 + 220 x + 26 x13 + 228 x + 27 x13 + 229 x + 28 x13 + 218 x + 29 x13 + 229 x + 210 x13 + 227 x + 211 x13 + 216 x + 212 x13 + 227 x + 213 x13 + 225 x + 214 x13 + 229 x + 215 x13 + 225 x + 216 x13 + 223 x + 217 x13 + 227 x + 218 x13 + 223 x + 219 x13 + 221 x + 220 x13 + 229 x + 221 x13 + 230 x + 222 x13 + 219 x + 223 x13 + 230 x + 224 x13 + 228 x + 225 x13 + 217 x + 226 x13 + 228 x + 227 x13 + 226 x + 228 x13 + 230 x + 229 x13 + 226 x + 230
Table 16: Irreducible trinomials of the form x13 + 2s1 x + 2s0 for 0 ≤ s1 , s0 ≤ 30 that can be used for construction of the finite field GF (q 13 ) where q = 231 − 1.
21
