CC BY
7
Дослиджено теоретичну складтсть виконання квазжвадратичних алгоритмiв обчислення кшько-стi точок елттичних кривих, визначеног над двш-ковим полем. Проведено експериментальний аналiз таких методiв з використанням програмног моделi. Описано необхiднiсть у збшьшенш розмiрiв загаль-носистемних параметрiв национального стандарту електронного цифрового тдпису. Наведено пришвид-шений метод за часом свого виконання, що може бути використаний для модертзацп стандарту
Ключовi слова: метод Сато, метод Харлi, порядок елттичних кривих двшкове поле, ^д Фробешуса □-□
Исследована теоретическая сложность выполнения квазиквадратичных алгоритмов вычисления количества точек эллиптической кривой, определенной над двоичным полем. Проведено экспериментальный анализ таких методов с использованием программной модели. Описана необходимость в увеличении размеров общесистемных параметров национального стандарта электронной цифровой подписи. Приведен ускоренный метод по времени своего исполнения, который может быть использован для модернизации стандарта
Ключевые слова: метод Сато, метод Харли, порядок эллиптической кривой, двоичное поле, след Фробениуса
UDC 004.056
|DOI: 10.15587/1729-4061.2017.95194]
EXAMINATION AND IMPLEMENTATION OF THE FAST METHOD FOR COMPUTING THE ORDER OF ELLIPTIC
CURVE
I. Gorben ko
Doctor of Technical Sciences, Professor Department of Security of Information Systems and Technologies V. N. Karazin Kharkiv National University Svobody sq., 4, Kharkiv, Ukraine, 61022 Е-mail: GorbenkoI@iit.com.ua R. H a n z i a Postgraduate student Department of Information Security Technologies Kharkiv National University of Radio Electronics Nauky ave., 14, Kharkiv, Ukraine, 61166 Е-mail: roman.ganzya@gmail.com
1. Introduction
The quality of providing electronic trusted services, primarily in the functioning of electronic government, is impossible without the use of electronic digital signature (EDS). Its application makes it possible to provide users with such services as integrity, authenticity, irresistibility and authenticity of the source of information, by which it was created. In Ukraine, to solve the indicated tasks, national standards are used - DSTU 4145-2002 and DSTU 7564-2014. In addition to the national standards, such harmonized ones are also used as, for example, DSTU ISO/IEC 14888-3 and DSTU ISO/IEC 9796-3. As shown by the analysis, cryptographic strength of the specified standards largely depends on the properties of general system parameters that are employed, the possibilities of expanding their dimensions and operative generation and replacement. The aforementioned becomes especially relevant in the case of emergence of quantum computers and their application for conducting cryptanalytic attacks, for example, to solve the problems on full disclosure and protection from collisions [1].
Therefore, absolutely relevant is a problem task on the prompt generation and application, with the possibility of replacing the tuples, of strong system-wide parameters for EDS. The EDS used in Ukraine are implemented based on cryptographic transformations in the group of points on elliptic curves. At present, the specified task of rapid replacement of general settings and dimension of these parameters
is not solved and there is an essential need to consider and resolve it [2].
2. Literature review and problem statement
The national standard of EDS DSTU 4145-2002 was adopted in 2002 and contains system-wide parameters with size up to 431 bits. At that time, it was considered that such parameters would be able to provide the necessary level of the cryptographic transformations stability on elliptic curves (EC) for decades. However, as shown in studies [1], the current state of technology development in the direction of constructing a quantum computer and existence of the Shor's quantum algorithm [3] open up new threats for the existence of modern cryptography, especially asymmetric, including that on EC.
The National Institute of Standards and Technology in the USA (NIST) plans to adopt in the next three to five years the first standard of asymmetric crypto-transformation. However, it is more likely that before its application there will be a "transition period" when both existing and post-quantum EDS are used [2]. Nevertheless, to ensure cryptographic strength of the existing EDS, for example, of those specified in DSTU 4145-2002, it is necessary to increase the dimensions of the system-wide parameters. Such an increase will not significantly improve stability against quantum cryptanalysis. However, for running such an analysis, a true quantum computer is required with a large number
of qubits, which will not be created for many years to come. Therefore, the use of EDS in line with DSTU 4145-2002 with an increased order of the base point to 1000 and more bits can be a guarantee for the cryptographic stability for a significant number of years [4].
The main problem on forming the general systems parameters for use in the groups of points on EC comes down to a complicated, in terms of calculation, task - to compute the number of points on EC [4]. The problem to calculate the number of points on EC is a non-trivial one and, at present, in Ukraine there are no data in the available sources on the execution order and essence of this stage. However, articles [5-13] present an overview and proof of mathematical methods that can be employed to count the number of points on EC.
The algorithm for computing the trace of Frobenius endo-morphism, presented in paper [9], is the first p-adic algorithm that was proposed to resolve this type of problems. Earlier, to calculate the order of EC, they used l-adic algorithms [5]. The article [6] can be considered a basis in the direction of calculation the order of EC for binary fields. Author in paper [6] proposed a modification of the method so that it can be used not only for the characteristic 3 and larger. The computational complexity of the method described there is rather high. In particular, article [8] presents a detailed analysis of existing methods for the calculation of zeta-functions and author of the work proposes certain improvements to the algorithm from paper [9]. All the improvements to the method from article [9] are based on the reduction of its computational complexity. In paper [10], authors propose certain mathematical improvements that will be analyzed in detail later in the present study. Thus, it is proposed to use another basis, analytical method for normalization and the application of inverse Frobenius substitution, in contrast to articles [6, 8, 9], where the ordinary substitution was employed. In [12], author proposes a certain improvement to the algorithm based on arithmetic and geometric progression, which is based on using other modular polynomial. Further, in the present study, there will be a detailed analysis of articles [10, 12] and we shall experimentally test the effectiveness of such improvement. Papers [11, 13] describe a variety of solutions to the Artin-Schreier equation for various bases. In [11], for optimal polynomial basis, in [13], for the Gaussian normal basis. A detailed theoretical and practical analysis of these studies will be presented further in our research.
Paper [14] analyzed the international standard ISO/ IEC 15946-5, which defined the methods for generating EC. Research in [14] demonstrated that for the finite field F(2m), significant advantages by the criterion of time complexity (performance speed) are displayed by the methods that are implemented based on the p-adic numbers. In this case, mathematics in the ring of p-adic numbers needed for constructing a program model was described in articles [7, 8, 15]. As the main criterion in the present study, we shall use the criterion of time complexity (performance speed) in the calculation of the number of points on EC and the guarantee of properties of the appropriate system-wide parameters. We shall also subsequently compare efficiency of the received results with the fundamental results that are presented in papers [8, 14].
3. The aim and tasks of the study
The aim of present research is to determine and implement optimal (by performance speed) algorithm for the
computation of EC order, which can be applied to increase the base of system-wide parameters of the Ukrainian EDS standard.
To achieve the set aim, the following tasks were formulated and resolved:
- analysis and substantiation of selecting the most promising methods for computing the EC order out of those existing;
- selection by the criteria of minimization of the computational complexity in the algorithm for generating system-wide parameters of high and super-high levels of stability;
- development of a program model and conducting a program simulation, as well as a comparison of theoretical and obtained experimental results.
4. Materials and methods for examining the algorithms to compute the order of elliptic curves
4. 1. Overview of canonical lift of EC
Article [14] presents a general approach to count the number of EC points (EC order) E, assigned by equation:
y2 + xy = x3 + a6, with element a6 over Fq with j-invariant j(E) = a-1
and q=pn. Such a general approach includes phases of canonical lift of EC, conducting the normalization of coefficients of lifted curve to the base field and the computation of trace of Frobenius endomorphism. Further attention will be paid to the optimal methods of canonical lift of EC that possess quasi-quadratic complexity of their execution or close to such. The stage of normalization is described in more detail in papers [7, 8] and will be described in the subsequent studies on the given subject.
The first algorithm for canonical lift of EC using the p-adic methods was proposed in article [9] (described in detail with certain modifications in paper [8]) and is based on the fulfillment of the Lublin-Serre-Tate theorem.
Theorem. Assume E is not a supersingular elliptic curve over the field Fq. Then, up to isomorphism, there is one curve s, defined above Zq, such that:
1. EC equation s is equal to EC equation E by module 2.
2. End(e) s End(E).
A corollary to this theorem is the following commutative diagram:
M-1
* £n
i ^
E _
1...... y °0
i y i y iy
<Pp,<H
-> E...
-^CP
(2)
where si is the canonical lift E and 9pI is the corresponding lift ¿p,i.
Such isogeny of degree 2 allows us to associate modular equations also of the 2 degree j-invariants of canonically lifted curves:
®2(Ke),Z№))) = 0, (3)
where 02(X, Y) is the symmetric bivariate polynomial:
e
O2 (X, Y) = X3 + Y3 - X2 Y2 + 24 ■ 3 ■ 31 (X2 Y + XY2) -
-243453 (X2 + Y2) + 3453 ■ 4027XY +
+283756 ■(X + Y)- 2123959. (4)
A more detailed description of the canonical lift is given in articles [8, 9]. In the following chapters, we shall demonstrate the application of basic properties of canonical lift, as well as their evolution to rapid counting of the number of points on EC.
4. 2. The Satoh-Skjernaa-Taguchi method (SST)
In the beginning of counting the number of EC points, j(E) in Fq/F4 is known and it is required to compute j(s) according to (3). The value of j(s) can be calculated by receiving in turns one by one bit of the entire magnitude: assume that we know J as:
J = j( e)mod2k, then we can record computation:
j( e) = J + 2ke and insert it formally into equation:
O 2(j(e), X(j(e))) = 0,
O p(X,Y) = (Xp - Y)(X - Yp), it follows that:
!^(J,5:(J)) + 2k = 0,
e* * 0,
module 2. Finally, since Z(e)=e2mod2, we obtain:
, o 2(J,z(J))
(6»
2k f^o»»
mod 2.
(7»
Taking the unique root of power 2 (e£Fq) leads to a better approximation of j(s) that is assigned by:
J + 2ke = j( e)mod2k+1.
In paper [10], authors proposed a different approach to solve such problem (difficult in terms of computational complexity) as finding the root for the situation of degree 2.
To avoid having to find the root, it was proposed to replace the system of equations in (2) with such J=j(s) by module 2 with:
hence, we obtain equation e by module 2. Thus, receiving one bit when approximating to j(s) [8].
In a general case, this method can be used for different characteristics of finite fields, but in the present study we shall carry out the adaptation towards extending the finite fields of characteristic 2. This is due to the existence of acting standard on the electronic digital signature based on elliptic curves DSTU 4145-2002. In a general case, expressions of type:
®2(Z-1(J),J) = 0 and J = j(E)mod2.
(8»
In this case, expression in (7) can be converted into an other form and thus avoid having to calculate the root:
_ ®2(2-1(J),J)
2k
mod2.
(9»
It is necessary to note that:
j(£) = J + 2ke,
take the form:
j( e» = J + pke
and others on the analogy.
For a more accurate tuning, let us take its decomposition:
O 2(j(e), Z(j(e))) = 0,
in a Taylor series:
0 = O2(j(e), Z(j(e))) = O2(J + 2ke, Z(J + 2ke)) = = O2^J,Z^J))+ 2ke* ||(J,:£(J))+
+2ke* dOYr(J, 2(J)) + 22kr(e, 2(e)),
with such element r that reZq[X,Y]. In this equation:
O 2(j(e), X(j(e))) = 0,
(5»
by module 2 in degree k, and hence we can divide expression by 2k and obtain ratio e by module 2. In addition, from the Kronecker relationship:
O2(X-1(J),J) = 0mod2k, and we only need to compute the inversion:
dOY2(^-1(J),J),
by module 2. The only problem that remains is the calculation of inverse Frobenius substitution.
Authors of the national standard of electronic digital signature DSTU 4145-2002 propose to use the representation when f(x)£Zq[x] is used as a module. It means that such a polynomial is the discharged one, irreducible to Fp[x] of degree n. Using such a representation makes it possible to effectively perform algebraic operations in a ring, especially the operation of taking by the module. However, for the operation, which is necessary to compute (7) or (9), in particular calculation of Frobenius substitution value, it is very slow.
To reduce computational complexity when calculating the Frobenius substitution, authors the SST method in article [10] proposed to use different representation of the ring of p-adic numbers. Such representation makes it possible to significantly accelerate computing the Frobenius substitution, but requires certain precomputations.
The essence of using alternative representation is in the following. Assume 6 e Fq is zero in f(x)mod p, in other
words, Fq = Fp[9]. It is obvious that 6 is the (q-1)th roots of unity because 6eFq. Assume 0 is the Teichmuller lift from 6, in other words, the (q-1)th root of unity in Zq and 6 = 6 mod p, then we define m(x)£Zp[x] as a formative polynomial. In this case, m(x)=f(x) by modulo 2, and hence it follows that it is possible to assign a ring of p-adic numbers (Qq as Qq[x]/m(x)) [8]. More details on the procedure for calculating the Teichmuller module can be found in [7].
Such a representation of the ring allows us to reduce the computation of Frobenius substitution. Since Z(0)=0p by module 2 and Z(0) is the (q-1)th root of unity in Zq, the authors conclude that:
2(6) = 6p mod p.
Therefore, efficient computation of the Frobenius substitution is reduced to:
2. y = jmod2;
3. fori = 2 to N do
3. 1. x = 2-1(y)mod2i; (13)
3. 2. y = y - d* O2(x,y)mod2i;
4. return y.
Substitution of lift step of the j-invariant with algorithm (13) in the Satoh algorithm will yield the acceleration of computing the EC order by about 5 times. Authors of the SST method pointed that at every step of such algorithm (13), it is necessary to recalculate 02(x,y), though values of x and y at step i+1 are very close to the values that are used at step 1. Therefore, we can compute:
y = j(E)mod2W
and
I Ia,6' = Ia,6'p,
(10)
the result of the calculation in this case should always be reduced by module m(x). For a binary field, characteristic p takes value 2.
Returning to (9), in particular to the problem on computing the inverse Frobenius substitution and employing new presentation of the ring, such computation takes the following form:
apk+j9
Cj(6),
x = 2-1(y)mod2W for certain W. Assume m > 1 and i > 0:
O2(x + 2mW+' Ax,y + 2mW+' Ay,) = O2(x,y)+
+2mW+' (^ddxL(X'y)Ax + ddY2(X'y)Ay) m°d2(m+1)W, (14)
(11)
i n-1 A p-1 (
I-1 lI^ = I I
V i=0 ) j=0 V0<pk+j<n
where
Cj(0) = I-1(ej)=ejpn-1.
This value can be precomputed (as well as the value of the Teichmuller module). Thus, if we perform all the necessary precomputations, then finding the value of inverse Frobenius substitution, in terms of computational complexity, will be simple. Interpretation of formula (11) to the binary field will take the following form:
I-1 (E^ 1 = E a2k6k* X a2k+16k*C1(6). (12)
V i=0 / 0<2k<n 0<2k+1<n
For a binary field, actually, calculating the inverse Frobe-nius substitution does not require a large amount of computations, and precomputations come down to finding only one element:
C1(6) = 62n-1,
that is also a trivial task.
If we take an idea to use in particular the inverse Frobe-nius substitution and, accordingly, other space, and combine with the Satoh algorithm, then it is possible to represent the idea of the SST method ("grey" method) as follows:
Algorithm. Lift_J_Naive
Input: j-invariant j £ F2n\F22 and accuracy N
Output: J £ Zq, J=j mod 2 and &2(£-1(J), J) = mod 2N.
1. d = (d|Y2( 2-1(j),j)1 1mod2;
so it remains to find out the values:
/ x a°2 / X
—-(x.y) and —-(x.y). 9X V y' 9Yv y'
Thus, it is possible to use the indicated equation to compute:
O 2(x + 2mW+i Ax, y + 2mW+i Ay),
from value 02(x,y) while i<W. This idea is represented in the following algorithm, which conditionally consists of two parts: in the first, we compute y = j(E)mod2W using algorithm (13), and in the next, we use expression from (14). Such expression is used to incrementally update the value of 02(x,y) without its recalculation at each step [10].
Algorithm. LiftJ
Input: j-invariant j £ F2n\F22 and accuracy N Output: J £ Zq, J=j mod 2 and &2(£-10), J) = mod 2N.
1. d = 2-1(j),j)l mod2;
I dY
2. y = jmod2;
3. fori = 2 to W do
3. 1. x = X-1(y)mod2j;
3. 2. y = y - d* ©2(x,y)mod2';
4. x = Z-1(y)mod2W;
5. D = ^(x,y)mod2W;
x 9X v y'
6. D = ^2(x,y)mod2W;
y aY v 'y> '
7. form = 1 to[(N- 1)/Wjdo
7. 1. x = X-1(y)rnod2(m+1)W; 7. 2. V = 02(x,y)rnod2(m+1)W; 7. 3. for i = 0 to W-1 do 7. 3. 1. Dy = -d * 2-(mW+j) V mod 2;
7. 3. 2. Ax = Z-1(Ay) mod2W-i;
7. 3. 3. y = y + 2mW+' Aymod2(m+1)W;
7. 3. 4. V = V + 2mW+'(DxAx + DyAy)mod2(m+1)W;
8. return y.
Authors in paper [10] proved that the complexity of calculation in this method, algorithm in (15), for W--n^/a+H) and N-n/2 is O(n2^+1/(1+H)). In practice, authors also recommend using such a value for W that is multiple to the word size of central processor. Of course, the above-described complexity does not account for the time to perform precomputations for the polynomial, which assigns field m(x), also the value Cj(0) from expression (11). One has to point that the complexity of computing C1(6) = 62 is O(n2H+1), demonstrated in [8].
In general, the full version of this algorithm can be shown in algorithm (16). This version is the full version of the SST algorithm and includes the following steps:
1. Computing a polynomial that assigns the field (can be recalculated).
2. Computing element Cj(0) for finding the inverse Frobenius substitution (can be precomputated).
3. Lifting the j-invariant of curve to the required accuracy.
4. Finding value Co from expression:
Vp*(x) = c0T0 + O(T2),
This step for the characteristic of field 2 is described in more detail in articles [8, 10].
5. Normalization of coefficients of value from step 4 and finding a trace of Frobenius endomorphism.
6. Obtaining value of the EC order as #E(Fq)=1+q±t, where t is the value obtained in the previous step.
Algorithm. SST
Input: Elliptic curve E: y2 + xy = x3 + c over F2d
Output: Number of points on curve E(F2d)
' d _
1. N =
-13;
2. M = N -10;
3. j = j(E»;
4. m(x» = GenTeichmullerModule(f(x»»;
5. C = GenC1(x»mod2N;
6. J = Lift_J(j»mod 2N;
7. J'=Z-1(J»mod2N;
8.
Function of the Teichmuller module generation: m(x) = GenTeichmullerModule(f(x)), described in detail in [7, 10], and by function:
C = GenC1(x), we mean the computation: C1(6) = 62"-1,
to find the value for the inverse Frobenius substitution (12). Function Lift__J(j) in the 6th step of algorithm (16) is described in more detail in (15). The functions of finding a square root in the ring of p-adic numbers used in step 12, as well as other mathematical functions that are defined in the ring of p-adic numbers, are described in detail in [7]. The operation of normalization in the 13th step of the algorithm implies the recovery of trace of Frobenius endomorphism [15]:
t = H21(co» = NQq/Qp (co »modq.
(17»
A fast algorithm for computing the norm and proof of its work was proposed by authors of the SST method in the same paper [10]:
NQq/Qp (a» = eXP(TrQq/Qp (l0g(a»»».
(18»
Z = --
(J'2 + 195120J' + 4095J + 660960000» 8(J'2 + J'(563760 - 512J» + 372735J + 8981280000»
9. T' = ((12Z2 + Z)(J -1728) - 36)mod2M;
10. CN = (J - (504 + 12096Z)T)mod2M;
11. CD = (T(240T + J))mod2M;
12. c = Sqrt(CN/CD)mod2M-1;
13. t = Norm(c)mod2M-1;
14. if t2 > 2d+2 then t ^ t - 2n-1;
15. return 2d +1 -1.
The algorithm presented requires certain explanations and clarifications, in particular j(E) = 1/c for the curve equation:
E : y2 + xy = x3 + c over F2d.
This algorithm works fast for finite fields of characteristic 2, but for finite fields of other characteristics its computational complexity increases significantly. This is primarily caused by the complexity of calculating the logarithm in the ring of p-adic numbers. It should also be noted that this algorithm in the wording proposed in [10] for the binary fields will work only for the sparse generating polynomials. For the polynomials that form the Teichmuller basis, expression from (18) with the possibility of rapid computation of the logarithm will not work out. This causes the need to transform element c after step 12 of algorithm (16) back to the optimal polynomial basis.
4. 3. The modified Satoh-Skjernaa-Taguchi method (MSST)
Paper [12] proposed to combine the idea of arithmetic-geometric method and the Satoh-Skjer-naa-Taguchi method. The AGM method was an-mod2N; (16) alyzed in detail in article [14]. It is worth noting that for the MSST work, one should use a one shift variation of the AGM method, though in terms of computational complexity, one shift and two shift AGM are almost indistinguishable. However, using the one shift variation implies applying the AGM-sequence (ak ,bk )k=0, in which:
ak = bk = 1mod4, ak = bkmod8, in the next variant, Xk=ak/bk, which matches EC: E. : y2 = x(x-1»(x -Xk».
(19»
As each preceding AGM-sequence makes it possible to compute the next one, in other words, such sequence is iterative:
k+1
) = ((ak + bk)/2^>/iA),
then we may represent the iterative function of the one shift AGM-sequence in the form
2Jh ^ t+-| =^L~L [12]. k+1 1 + Xk L J
Initialization of the one shift AGM-sequence is performed as follows:
^ = (1 + 8c)mod16,
(20)
E(X,Y) = Y2(1 + X)2 - 4X = 0.
(22)
The above representation allows us to use the condition of modular polynomials (Lublin-Serre-Tate theorem) in the form:
E2(X, Z(X)) = 0mod2k+3.
It should be noted that both partial derivatives are equal to zero by module 2 in this expression, which is why it is not possible to directly use the SST algorithm. To eliminate such a problem, [12] proposes the following substitutions:
X ^ 1 + 8X,
Y ^ 1 + 8Y
and, as a result, we shall obtain a modified modular polynomial for finite fields of characteristic 2 in the form of:
equal to unity that satisfies the requirements of the SST algorithm [10].
The last thing that is necessary to do is to align the expression from which a Frobenius trace is computed, for the two-shift AGM sequence it takes the following form:
TrF = tk + q/tkmodk
(25)
from
where c = c mod 2 is the free EC coefficient. Gaudri also proves that sequence Xk+i converges similar to ak/bk, it is implied that:
Xk+1 =X(Xk)mod2k+3,
and if we substitute this value into the value of iterative function of the AGM-sequence, we shall obtain:
E(Xk)2(1 + Xk)2-4Xk = 0 mod 2k+3. (21)
The above-indicated expression is solved with the use of the main idea of the SST algorithm. Paper [12] presents the solution to this problem through the change of expression for modular polynomial that is used for EC lift (described in (3)).
Assume E(X,Y) is the module expression of AGM, then:
tk = NQq/Qp (ak / ak+1 )
substituting in this expression the values that were accepted for the one-shift AGM, in other words:
^k = ak/bk, ak+1 = (ak+bk)/2
and
= 1 + 8c,
we shall obtain:
tk = nq,/qp
1
1 + 4L
(26)
Upon presenting the SST method, the algorithm described above does not need special explanations. If one compares two similar methods of SST and MSST in terms of computational complexity, the MSST method is usually works faster. For W=n^/(1+^) and N=n/2, its complexity is O(n2^+0,5). This is caused by the fact that the modular polynomial used in the MSST method requires 1 multiplication and 1 squared raising. The SST method needs 3 multiplications and 2 squared raisings (meaning in the ring of p-adic numbers). From the point of view of the spatial complexity, the two methods require similar resources. Spatial complexity for them is O(n2). Although the constants that are used in the MSST modular polynomial are significantly lower than those in the SST method. The phase of precomputations for two algorithms is the same and requires O(n2^+1) computing resources for its implementation [8].
4. 4. The Harley method
Article [13] proposed a p-adic algorithm for finding the EC order without precomputations and with computing complexity O(nMogn). The main idea of this method was applying the method for solving equations by Artin-Schreier, that is, equation of the form:
E = (X + 2Y + 8XY)2 + Y + 4XY = 0,
(23)
that can be solved when X is known and it is not equal to zero by module two, so that Y = Z(X). In this case, the partial derivatives are equal to:
—(X,Y) = 2(X + 2Y + 8XY)(1+8Y) + 4Y, (24)
9X
dE
—(X,Y) = (4(X + 2Y + 8XY) +1)(1 + 4X), 9Y
it is not difficult to confirm that each partial derivative by X is indeed equal to zero by module two, and by Y is
xp-x + a = 0,
with element a£Fq (by the field characteristic - p). Positive version of Hilbert space claims that such equation has the solution in Fq only in the case when:
TrFq/Fp (a) = ^
Since: X(x) = xp,
then the given equation for Zq takes the form:
E(x) - ax - b = 0, (27)
with a, beZq. Let us define:
E k(x) = akx + bk, for all k=2,...,n. As E n(x) = x,
the authors conclude that the above-given equation is assigned as:
bn/(1 - a„),
and for computing ak, bk, the authors propose the following formula
Ek+'(x) = E'(akx + bk) = E'(ak)(a,x + b.) + E'(bk). (28)
It is the above-given expression that is used to solve the Artin-Schreier equation. Next, to lift EC, a generalized Newton's algorithm is used [13], the very basic idea of lifting remains the same as in the SST algorithm (2). However, [13] proposed these parameters (solving the Artin-Schreier equation and the generalized Newton's algorithm) for the Gaussian normal basis only.
[11] presents a similar variant for solving the Artin-Schreier equation and a different generalized Newton's algorithm that can be used for the optimal polynomial basis and the Teichmuller basis. Assume we have the following problem: given
O(X,Y) eZq[X,Y],
and it is required to find the root x£Zq such as:
O(X,E(X)) = 0. (29)
Assume we know that:
xm = xmodpm
and let it be
8 m = (x - xm)/pm,
then the decomposition in a Taylor series around xm for p=2 yields:
0 = O2(x, E(x)) = O2(xm + 2m 8m, E(xm + 2m 8m)) = = O2(xm, E(xm)) +2m (8m Ax + E(8m)Ay) mod22m, (30)
- p2(x,E(x)) 2m
Assume:
= 8Ax + E(8)Aymod2m.
(31)
with
Ax = d°2(xm, E(xm))mod2m 9X
and
Ay = !y (xm,E(xm))mod2m.
Hence, it follows that 8m has a solution of the following type:
k = ord2(Ay),
and then if ord2(Ax) > k and
ord2 (xm, E(xm)) > k + m
and m>k, we obtain the following equation:
aE(8) + p8 + g = 0mod2m-k, (32)
coefficients a, p, y£Zq and a are the unity in Zq. As the Frobenius substitution maintains the values, then the above-given equation computes Sm mod 2m-k unequivocally and, in addition:
x = xm + 2m 8mmod22m-k.
We shall assume that there is an algorithm that returns the null value 5t' from equation (32) with accuracy t' = [ t / 2] (half precision), then there is a possibility to employ the same algorithm for computing St with accuracy t (full precision). Let us substitute:
8t = 8t. + 2'At
we shall receive in (31):
aE(At)+bAt + aE(8tHp8 t+g = 0mod2"', (33)
that is why, since t -1' < t', then it is possible to apply the same algorithm for computing At by module 2t-t' and, therefore, we may receive the value St. Thus, it gives us a recursive algorithm that allows us to find value (32). If we assume that:
ord2(Ax) > ord2(Ay),
then it follows that ord2(P) > 0 and everything is reduced to solving:
aE(8) + g = 0mod2.
And since a is the unity, it allows unambiguously computing 5 mod 2. To calculate the Artin-Schreier equation, it is proposed to employ the following algorithm [11].
Algorithm. Artin-Schreier-root
Input: a, p, yEZq, aEZxq, ord2((¡)>0, accuracy N.
Output: xEZq so as aZ(x)+Px+y=0mod 2N.
1. if(N == 1) then
1.1. x = (-g / a)1/2 mod 2;
2. else
2. 1. N' = [N/2];
2. 2. M = N - N'; (34)
2. 3. x' = Artin - Shreier - root(a,P, g,N'); 2. 4. g' = (aE(x') + Px'+g)/2Nmod2M; 2. 5. A' = Artin - Shreier - root(a,P, g',M);
2. 6. x = x'+ 2n'A 'mod2N;
3. return x.
Paper [8] demonstrates certain optimization regarding computing in step 1. 1. of algorithm 4 (34), in other words, root calculation takes the following form:
p-1
/ n-1 \ 1/p E^ =E E aj
j=0 V0<pk+j<n
V i=0 /
where
Cj(6),
(35)
Cj(6) = (6j)1/p = 6jpn-1 and version for the binary field takes the form:
( n-1 \ 1/2
I^6i = E a2k6k* X a2k+16k*C1(6),
V i=0 1 0<2k<n 0<2k+1<n
(36)
where
C1(6) = 62-.
Assessment of the complexity of the above-given algorithm is also provided in article [11] and amounts to:
O((nN)MlogN).
The complexity of the algorithm as a whole is based on recursive challenges in steps 2. 3. and 2. 5. of algorithm 5 (34), multiplication and computation of the Frobenius substitution in step 2.4. The complexity of computing the Frobenius substitution is O((nN)^) bit operations [8].
The above-given algorithm allows us to effectively solve equation (33), and to solve the main problem outlined in (30) and reduced to (32) in paper [11], proposes to use a generalized Newton's lift somewhat different from the one described in article [13].
Algorithm. Generalized-Newton-Lift
Input: Modular polynomial x0 £Zq, which satisfies
O2(x0,E(x0))=0 mod 2k+1 and f^0-(x0,E(x0))|>k, k= \ V dX 1
=ord.
a®,
x0, E(x0)) I, accuracy N.
9Y
Output: xn£Zq, 02(xn,Z(xn))=0 mod 2N+k and xn= =x0 mod 2k+1.
1. if(N < 2k +1) then
1. 1. x = x0;
2. else
2. 1. N' = [" N/2] + k;
2. 2. M = N - k; (37)
2. 3.
Ax = dOL (x',y') mod2N) / 2kmod2M;
x' = Generalized - Newton - Lift(x0 ,N');
2. 4. g' = E(x')mod2N;
2. 5. V = (O 2(x', g ')mod2N)/2N'mod2M;
2. 6. V = (O2(x', g ')mod2N)/2N'mod2M;
2. 7. Ay = dOL (x',y') mod2N) / 2kmod2M;
2. 8. A' = Artin - Schreier - Root( Ay, Ax, V,M);
2. 9. x = x'+ 2M A 'mod2N;
3. return x.
The complexity of algorithm (37) is similar to algorithm (35). The main complexity is based on the recursive challenges of the algorithm from (35).
Multiplication of two integers that consist of n bits is performed in O(n^) operations, where ^ is the constant that defines the period of multiplying two m bit integers with time complexity O(m^). Thus, for classical algorithms of multiplication, values |j=2, and for the fast Karatsuba algorithm, |j=log23.
It should be noted that each of the presented methods finds only the order of the curve and does not tackle the issue on the possibility of its use in cryptographic transformations. Based on this, upon computing the order of the curve, it is necessary to verify the feasibility of its application in the cryptographic systems. And to choose the criteria that make it possible to select elliptic curves for constructing system-wide parameters at the required level of stability.
5. Results of exploring time complexity of the methods for computing the order of elliptic curves
When exploring the mechanisms for constructing strong cryptographic system-wide parameters of EC, an important criterion is the time required for the construction of such parameters and the dimensionality of field, over which the curve is defined. Parameters obtained in this way will be suitable for their further use in the EDS mechanisms according to DSTU 4145-2002 (or similar, in other words those that employ a binary field). The Ukrainian standard sets the parameters of size up to 431 bits while the standard FIPS 186-3 contain system-wide parameters up to 521 bits. Important is the estimation of time needed to build strong cryptographic parameters at the super high level of stability (509 < #E < 1031, bits) [4].
The largest complexity and resource consumption when generating the parameters is displayed exactly by the step in counting the points on EC, the theoretical information about this stage is presented in the previous chapters of present article. Those chapters also describe theoretical evaluation of complexity in the phase of EC lift for all presented algorithms. The complexity of computing the norm is given for only one of the existing algorithms that was used to calculate the order of EC.
For the canonical lift, we applied the SST method [10], the MSST method [12] and the Harley method [11] (described in chapters 4.2-4.4 of the present article). For the normalization, an analytical method (proposed in paper [10]) was employed.
To count the number of points on EC, we developed a software tool in the C++ language using the library NTL and gmp. A research into algorithm execution time was carried out in the program that was compiled using gcc 4.84 in the operating system Ubuntu 14.04 (USA) and the processor CPU Intel Core i5-2300 (USA). As all the operations for this class of algorithms are conducted sequentially, the par-allelization of implementation of the algorithm is impossible, and the number of cores in the processor will not change execution time of the algorithm.
Table 1 gives the period of computing the lift phase for the SST, MSST, Harley methods and the norm. An analysis of Table 1 allows us to argue that there is an advantage of the MSST method over its standard modification in the SST method; the practical results of present study confirm theoretical research into the given methods. The Harley method is the most optimal among the three presented methods from the point of view of computational complexity.
Table 1
Computational complexity of p-adic methods
The degree of extension of the field d, bits Lift phase Normaliza tion by SST, s
SST, s MSST, s The Harley method, s
7 0,006265 0,003216 0,000925 0,000248
23 0,023238 0,010162 0,004397 0,00088
79 0,198349 0,077477 0,054023 0,007455
107 0,251051 0,123306 0,07934 0,010746
173 0,680349 0,400715 0,219232 0,028311
199 0,822255 0,501402 0,278724 0,035182
257 0,950504 0,698594 0,290267 0,054191
307 1,55134 1,17776 0,431548 0,088748
383 2,39215 1,91847 0,614475 0,133243
433 2,78954 2,245 0,766383 0,159998
503 3,97648 3,2511 1,00376 0,198481
601 7,15564 6,00326 1,60194 0,424093
709 10,44 9,02212 2,16152 0,532264
787 13,2826 11,6431 2,61117 0,649224
827 15,0722 13,1966 2,86254 0,715296
929 19,5872 17,5986 3,54738 0,920241
1021 24,9435 22,4928 4,35925 1,06235
1049 34,2923 30,0397 5,27848 1,7509
It can be argued that the canonical lift of EC using the Harley method is about six times more efficient than the SST method (the idea of which was employed in the Harley method). If we refer to the results of studying the AGM method and Satoh method in the previous article [14], then it can be argued that the Harley method is about 25 times more effective than the AGM method and is 75 times more efficient than the best modification of the Satoh method.
Fig. 1 shows a chart of dependence of the size of the field and the execution time for different methods of the canonical lift. A difference in the execution time between the algorithms of lift is observed for all the examined sizes of field extension.
......SST--MSST -The Harley method
40 -
35 -7
a, 30 -¿7
Td :
-g 20
r- — — r- — r- o> — — — m a> r- o> r- a> — r- cr, n-( r- r- a> r- — HI \Q O; rn '■€> 0> OO rn r- — in O ^ a> r- rr, r- c ] r- C] [-- ' ]
— — — c-] ci m m ^t >n '■o '■o r- r- oo oo a> a. o
Field size, bits
Fig. 1. Chart of dependence of computational period of canonical lift on the size of field extension for different algorithms
Fig. 2 shows a chart of dependence of the field size and the norm computation period for different methods of canonical lift of elliptic curves, although the same method of normalization was actually used.
......SST normalization — —MSST normalization
-Harley normalization
<>. so o m c o -t oc m r — 'r o. 'n o -t i m i c I [ c i r r i ^ ^ ^ r l C-l m n 't ^t 't ^ >0 h h B » O! Ö! o
Field size, bits
Fig. 2. Chart of dependence of the norm computation period for different size of the field extension
The chart in Fig. 2 demonstrates that such a method for computing the norm is efficient and results for all methods were about the same each time with a difference of 0.01 s.
6. Discussion of results of exploring complexity in the methods for computing the order of elliptic curves
The results obtained can be compared to the results received in article [8]. To count the number of points on EC, the author used the processor AMD XP 1700+ (USA) and the operating system Linux Redhat 7.1 (USA). The algorithms were written in the programming language C, and basic mathematical operations in the ring of p-adic numbers were performed on Assembler. Results of comparing the time indicators from [8] to the ones received in the present study are given in Table 2 (columns with data from article [8] contain the designation "B").
Table 2
Comparison of the time complexity of computing the order of curve
Extension of the field d, bits Total time for computing the order of curve
SST, s SST "B", s MSST, s MSST "B", s The Harley method, s The Har-ley method "B", s
144 0,44 0,13 0,24 0,06 0,17 0,06
168 0,63 0,26 0,37 0,08 0,2 0,08
192 0,74 0,29 0,45 0,13 0,26 0,12
240 0,77 0,65 0,55 0,28 0,23 0,25
288 1,3 0,72 0,97 0,39 0,37 0,38
336 1,83 1,17 1,43 0,64 0,49 0,6
384 2,4 1,76 1,91 0,97 0,6 0,92
480 3,44 3,56 2,87 2,03 0,91 1,87
Results in Table 2 demonstrate that the basic mathematical operations, written in the low-level programming language (implemented in paper [8]), yield large effectiveness when computing the order of EC. The most effective method for computing the order of EC, as shown in Table 2,
is the Harley method, in other words, research results of the present work and of the studies in article [8] regarding the method that is optimal in terms of computational complexity, coincide.
One should consider the size of 1031 bits, because this is exactly the size of elliptic curves required for the super high level of stability (512 bits for symmetric cipher). The total time for counting the number of points on EC of size 1031 bits for the Harley method at computing the norm by SST is approximately 10 s. This period slightly exceeds that of Table 1 because there are certain operations performed, which were not dealt with in the present work in detail. For example, present article does not address the process of polynomial generation for the field, converting elements before the normalization, etc.
It can be argued based on the presented findings that there is a convergence between the theoretical computational complexity of the examined algorithms and the compared experimental results. Similar to article [8], the present study observes reduction in the computational complexity from the SST, MSST methods to the Harley method. The obtained experimental assessments confirm the analytical complexity that was described by authors of the estimated methods in articles [10-13]. We can state that a combination of the Harley algorithm for EC lift and the normalization method from the SST algorithm are the best candidates to modify the Ukrainian EDS standard. Here by modification we mean an expansion in the base of general parameters. And, of course, the given combination is the best one only by the criterion of computational complexity.
The given article presents an optimal algorithm to count the number of points on EC that applies the Harley method [11] and the method for normalization form paper [10], as well as certain adaptation for a binary field used in the Ukrainian and world standards. Results of research into these algorithms demonstrated that they might be employed to modify the Ukrainian standard DSTU 4145-2002 in the direction of extending the number of general parameters and their size. The software model makes it possible to generate system-wide parameters at the super high level of stability in seconds.
In addition, present study shows the time required to generate parameters at the good level of stability (for example, 257 bits) - it is less than 0.5 s. At such characteristics, users of information systems can generate common parameters all by themselves prior to the phase of assigning keys between the parties. Previously, such a situation was impossible due to the high complexity of computations, but modern computing power and efficient mathematics provide the users with such possibility.
7. Conclusions
1. The conducted analysis of promising methods for computing the EC order revealed that to solve this problem, at present, the most efficient (in terms of computational complexity) is the Harley method. This conclusion was made based on the performed theoretical and experimental studies and comparison between the Satoh, AGM, SST, MSST and Harley methods.
2. Analytical complexity in the execution of the SST, MSST and Harley methods was demonstrated and, according to it, we determined the fastest (by execution time) algorithm for computing the order of the curve. It is demonstrated that the Harley method is the fastest, due to applying the Artin-Schreier equation for canonical lift of EC. By reducing the number of computations for the large size accuracy, the Harley method is more efficient than the SST and MSST methods. Present study shows that using the Harley method in practice makes it possible to accelerate the computation of EC order by approximately 7 times compared with the SST method.
3. Based on the explored data, we constructed a program model for the methods of canonical lift of EC and normalization. Development of a software model made it possible to perform experimental analysis of the examined algorithms. By the data obtained, the present work experimentally confirmed the quasi quadratic dependence of the field size, over which EC is defined, and the time required for the EC canonical lift.
References
1. Horbenko, Yu. I. Analysis of the possibility of quantum computers and quantum computings for cryptanalysis of modern cryptosystems [Text] / Yu. I. Horbenko, R. S. Hanzia // Eastern-European Journal of Enterprise Technologies. - 2014. - Vol. 1, Issue 9 (67). - P. 8-16. - Available at: http://journals.uran.ua/eejet/article/view/19897/18759
2. Hanzia, R. S. Analiz shlyakhiv rozvytku kryptohrafiyi pislya poyavy kvantovykh kompyuteriv [Text] / R. S. Hanzia, Yu. I. Horbenko // Visnyk Natsional'noho universytetu "L'vivs'ka politekhnika": Kompyuterni systemy ta merezhi. - 2014. - Issue 806. -P. 40-48.
3. Shor, P. W. Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer [Text] / P. W. Shor // SIAM Journal on Computing. - 1997. - Vol. 26, Issue 5. - P. 1484-1509. doi: 10.1137/s0097539795293172
4. Horbenko, I. D. Prykladna kryptolohiya [Text]: monohrafiya / I. D. Horbenko, Yu. I. Horbenko; KhNURE. - Kharkiv: Fort, 2012. - 868 p.
5. Schoof, R. Counting points on an elliptic curve over finite fields [Text] / R. Schoof // Proc. Journees Arithmetiques. - 1995. -Issue 93. - P. 219-252.
6. Skjernaa, B. Satoh's algorithm in characteristic 2 [Text] / B. Skjernaa // Mathematics of Computation. - 2003. - Vol. 72, Issue 241. - P. 477-488. doi: 10.1090/s0025-5718-02-01434-5
7. Handbook of Elliptic and Hyperelliptic Curve Cryptography [Text] / H. Cohen, G. Frey, R. Avanzi, C. Doche, T. Lange, K. Nguyen, F. Vercauteren (Eds.). - NW.: Chapman & Hall/CRC, 2005. - 807 p. doi: 10.1201/9781420034981
8. Vercauteren, F. Computing zeta functions of curves over finite fields [Text]: diss. for the degree of PhD / F. Vercauteren. - Katho-lieke Universiteit Leuven, 2003. - 195 p.
9. Satoh, T. The Canonical Lift of an Ordinary Elliptic Curve over a Finite Field and its Point Counting [Text] / T. Satoh // J. Ra-manujan Math. Soc. - 2000. - Vol. 15, Issue 4. - P. 247-270.
10. Satoh, T. Fast computation of canonical lifts of elliptic curves and its application to point counting [Text] / T. Satoh, B. Skjernaa, Y. Taguchi // Finite Fields and Their Applications. - 2003. - Vol. 9, Issue 1. - P. 89-101. doi: 10.1016/s1071-5797(02)00013-8
11. Harley, R. Asymptotically optimal p-adic point-counting [Electronic resource] / R. Harley // E-mail to NMBRTHRY list. -2002. - Available at: https://listserv.nodak.edu/cgi-bin/wa.exe?A2=ind0212&L=NMBRTHRY&F=&S=&P=7824
12. Gaudry, P. A Comparison and a Combination of SST and AGM Algorithms for Counting Points of Elliptic Curves in Characteristic 2 [Text] / P. Gaudry // Lecture Notes in Computer Science. - 2002. - P. 311-327. doi: 10.1007/3-540-36178-2_20
13. Lercier, R. Counting Points on Elliptic Curves over Finite Fields of Small Characteristic in Quasi Quadratic Time [Text] / R. Ler-cier, D. Lubicz // Lecture Notes In Computer Science. - 2003. - P. 360-373. doi: 10.1007/3-540-39200-9_22
14. Hanzia, R. S. Otsinka obchyslyuval'noyi skladnosti metodiv pidrakhunku kil'kosti tochok na eliptychniy kryviy [Text] / R. S. Hanzia // Systemy obrobky informatsiyi. - 2016. - Issue 8. - P. 92-99.
15. Satoh, T. Asymptotically fast algorithm for computing the Frobenius substitution and norms over unramied extension of p-adic number fields [Text] / T. Satoh. - Department of Mathematics, Faculty of Science, Saitame University, 2001. - P. 1-21.
Взв'язку з незадовшьною стштстю стандарт-них криптографiчних алгоритмiв з видкритим ключем до методiв квантового криптоаналiзу, проведено дослидження можливостi використан-ня постквантових криптографiчних алгоритмiв. Проведено порiвняльну оцтку таких алгоритмiв в залежностi вгд умов використання та проаналi-зовано переваги рiзних механiзмiв криптографiч-них перетворень, що е стшкими до методiв квантового криптоаналiзу
Ключовi слова: постквантовi криптографiчнi алгоритми, порiвняльна оцтка криптоалгорит-
мiв, критери порiвняння криптоалгоритмiв □-□
В связи с неудовлетворительной стойкостью стандартных криптографических алгоритмов с открытым ключом к методам квантового криптоанализа, проведено исследование возможности использования постквантовых криптографических алгоритмов. Проведена сравнительная оценка таких алгоритмов в зависимости от условий применения и выполнен анализ преимуществ разных механизмов криптографических преобразований, стойких к методам квантового криптоанализа
Ключевые слова: постквантовые криптографические алгоритмы, сравнительная оценка криптоалгоритмов, критерии сравнения криптоалгоритмов
UDC 004.056.55
|DOI: 10.15587/1729-4061.2017.96321]
EXAMINING A POSSIBILITY TO USE AND THE BENEFITS OF POST-QUANTUM ALGORITHMS DEPENDENT ON THE CONDITIONS OF THEIR APPLICATION
I. Gorben ko
Doctor of Technical Sciences, Professor* Е-mail: GorbenkoI@iit.com.ua V. Ponomar Postgraduate student* Е-mail: Laedaa@gmail.com *Department of Security of Information Systems and Technologies V. N. Karazin Kharkiv National University Svobody sq., 4, Kharkiv, Ukraine, 61022
1. Introduction
Due to the development of technologies for quantum computing and the introduction of quantum computer, there is a threat to the current state of protection of cryptographic systems with a public key [1]. With an advent of quantum computer that would have the volume of register required for the methods of quantum cryptanalysis, the stability of existing crypto algorithms will significantly degrade [2, 3]. This necessitates the creation of algorithms resistant to the methods of quantum cryptanalysis. The European project "New European Schemes for Signatures, Integrity, and Encryptions " (NESSIE) and the National Institute of Standards and Technologies (NIST) of the USA announced a start of
recruiting the applicants for the contest of post-quantum algorithms whose standards are planned to be adopted over 2020-2022 [4, 5].
A peculiarity of this task is that the contest will accept the algorithms whose cryptographic transformations are based on the latest information or insufficiently tested mathematical methods that will require considerable time to prove their stability in terms of quantum cryptanalysis. That is why the choice of the new standard will affect not only the algorithm that will be employed but also further development of the post-quantum cryptography.
Another feature is that the universal algorithms are lacking that can be used both for electronic signature (ES) and the encryption. Therefore, it is necessary for each of the security
