DOI QR코드

DOI QR Code

On the Use of Odd-Degree Isogenies for Implementing SIDH-Based Cryptography

SIDH 기반 암호 구현에 대한 홀수 차수 아이소제니 적용

  • Received : 2020.12.21
  • Accepted : 2021.01.25
  • Published : 2021.02.28

Abstract

In this paper, when SIDH is instantiated using only 3- and 5-isogeny, we demonstrate which curve is more efficient among the Montgomery, Edwards, and Huff curves. To this end, we present the computational cost of the building blocks of SIDH on Montgomery, Edwards, and Huff curves. We also present the prime we used and parameter settings for implementation. The result of our work shows that the performance of SIDH on Montgomery and Huff curves is almost the same and they are 0.8% faster than Edwards curves. With the possibility of using isogeny of degree other than 3 and 4, the performance of 5-isogeny became even more essential. In this regard, this paper can provide guidelines on the selection of the form of elliptic curves for implementation.

본 논문에서는 3차, 5차 아이소제니만 이용해 SIDH를 구현할 경우 몽고메리, 에드워드, 허프 곡선 중 어느 곡선에서 더 효율적인지 분석한다. 본 논문에서는 각 타원곡선의 형태에 대해 SIDH 암호를 구성하는 단위연산에 대한 연산량을 비교한 뒤, 홀수 차수만 활용해 SIDH를 구현하기 위해 소수와 파라미터를 설정하는 방법에 관해 설명한다. 본 논문의 결과 몽고메리와 허프 곡선에서 연산량은 유사하며, 에드워드 곡선보다 0.8% 효율적임을 알 수 있다. SIDH 기반 암호에 대한 다양한 파라미터 사용 가능성으로 인해 5차 아이소제니 구현은 필수적이므로, 본 논문은 이러한 SIDH 기반 암호에 대해 어느 타원곡선을 선택해야 하는지에 대해 가이드라인을 제공할 수 있다.

Keywords

I. 서론

대부분 공개키 암호의 안전성은 수학적 난제에 기반을 두고 있다. 현재 사용하는 RSA, DSA, ECC의 경우 각각 인수분해의 어려움, 이산대수의 어려움, 타원곡선 이산대수 문제에 기반을 두고 있다. 일반적인 컴퓨팅 환경에서 해당 문제들은 하지수시간 및 지수시간의 복잡도를 가지고 있어 안전하나, 만약 양자컴퓨터가 개발되어 Shor 알고리즘이 구현될 경우 해당 문제는 다항시간의 복잡도로 감소하여 더 이상 공개키가 안전하지 않게 된다 [16]. 한편, Shor 알고리즘이 구현 가능한 양자 컴퓨터 개발이 가시화되면서 이에 대응하기 위한 후 양자 암호 (Post-quantum cryptography, PQC)에 관한 연구가 활발히 이루어지고 있다. PQC란 일반 컴퓨터에서 구현할 수 있지만, 양자 컴퓨팅 환경에서도 안전한 암호를 의미하며, 기반문제에 따라 다변수 기반, 래티스 기반, 코드 기반, 해시기반, 아이소제니 기반의 5가지로 나눌 수 있다. NIST는 2016년에 PQC 암호에 관한 표준화 공모를 시작한다는 발표를 하였으며, 현재는 Round 3에 진입하였다. Round 3은 7개의 최종 후보와, 8개의 예비 후보 (alternate candidate)로 구성되어있다.

한편, 아이소제니 기반 암호는 2006년 Couveignes[3] 와 Stolbunov에 의해 처음 제안되었으며, 해당 기법은 오늘날 CRS라 칭한다. 하지만 제안한 암호에 대한 하지수시간 공격이 존재할 뿐만 아니라, 현실적으로 사용하기에는 속도가 매우 느려서 더 이상 연구되지 않았다. 2011년 De Feo와 Jao 에 의해서 Supersingular Isogeny Diffie-Hellman (SIDH)가 소개되면서, 아이소제니 기반 암호는 다시 주목을 받기 시작하였다 [10]. SIDH는 supersingular 타원곡선의 endomorphism ring이 가환이 아니라는 성질을 사용하여서, CRS 스킴 공격에 사용되었던 공격 기법이 적용되지 않아 가장 효율적인 양자 공격 기법이 지수시간 복잡도를 가지게 되어 안전하다. 또한, 기존 CRS 기반 암호보다 효율적이고 다른 PQC 암호들보다 키 사이즈가 작다는 점으로 SIDH 기반으로 설계된 key encapsulation mechanism 인 SIKE는 현재 Round 3의 예비 후보로 선정되었다 [1].

아이소제니 기반 암호의 구현은 크게 아이소제니 함수 구현과 타원곡선 연산 구현으로 나뉠 수 있다. 아이소제니 함수 구현 측면에서 바라보면, 해당 암호에서 사용하는 아이소제니 차수는 암호가 구현되는 유한체와 연관이 있다. SIDH 기반 암호의 경우 \(p=\ell_{A}^{e_{A}} \ell_{B}^{e_{B}} f \pm 1\) 형태의 소수를 사용하고, ℓA , ℓB가 사용되는 아이소제니 차수에 해당한다. 타원곡선 연산 측면에서 살펴보자면, 아이소제니 기반 암호의 기존 타원곡선 암호와 달리 고정된 하나의 타원곡선을 사용하지 않기 때문에, 타원곡선 연산을 타원곡선의 특성에 맞춰 최적화할 수 없다. 따라서 임의의 타원곡선에서 가장 효율적인 타원곡선 연산을 제공해주는 타원곡선의 형태를 사용해야 한다. 종합해보자면, 아이소제니 기반 암호를 구현할 때에는, 해당 알고리즘에 사용되는 아이소제니 차수의 구현이 효율적이면서 타원곡선 연산을 효율적으로 제공해주는 타원곡선의 형태를 사용해야 한다. 현재는 몽고메리 타원곡선을 이용해 대부분 아이소제니 기반 암호가 구현되어있다. 또한, SIDH 에 대한 state-of-the-art 구현도 몽고메리 타원곡선을 이용하고 있다.

하지만 최근에는 몽고메리 타원곡선보다 효율적인 다른 타원곡선의 형태가 존재하지 않을까에 관한 연구가 활발히 진행되고 있다. 첫 번째 후보로는 에드워드 곡선으로, 에드워드 곡선과 몽고메리 곡선이 birationally equivalent 하다는 점과, 두 타원곡선 사이의 변환이 효율적이라는 점으로 많은 연구가 진행되었다. 에드워드 곡선을 아이소제니 기반 암호의 사용은 Meyer 등에 의해 처음으로 제안되었다 [13]. Meyer 등은 [13]에서 SIDH 구현 시 아이소제니 연산은 몽고메리 곡선으로, 타원곡선 연산은 에드워드 곡선으로 구현하는 방법을 제안하였다. 이후, [12]에서는 아이소제니 연산이 에드워드 곡선에서 더 효율적이라는 점에 착안해, 타원곡선 연산은 몽고메리 곡선에서, 아이소제니 연산은 에드워드 곡선에서 수행하는 방법을 제안하였다. 특히 [11]에서 제안된 에드워드에서 홀수 차수 연산 공식이 몽고메리와 유사할 뿐만 아니라, 이미지 곡선의 계수를 복원하는 연산이 몽고메리보다 더 효율적이어서 에드워드 곡선은 더 주목받고 있다.

한편, 최근에는 허프 곡선을 이용한 연구가 활발히 진행되고 있다. 허프 곡선에서의 아이소제니는 [15]에서 처음으로 제안되었다. 하지만 타원곡선 연산과 아이소제니 연산이 비효율적인 관계로 주목받지 못하다가 최근 [6]와 [8]에서 효율적인 타원곡선 및 아이소제니 연산을 위한 좌표계가 제안되면서 다시 주목받고 있다. 특히, 해당 좌표계를 이용하면 3차와 4차 아이소제니 공식의 연산량이 몽고메리 곡선과 같아서 허프 곡선을 이용해 효율적인 SIDH 구현이 가능하다는 것을 보여준다 [8].

본 논문에서는 홀수 차수 아이소제니만 이용하여 SIDH를 몽고메리, 에드워드, 허프 곡선에서 구현한 결과를 제시한다. 일반적으로 SIDH는 3차 4차 아이소제니만 이용하는 파라미터만 제시되었다. 하지만 아이소제니 기반 그룹 키 합의를 구성할 때 [2] 사용되는 소수의 형태는 \(p=\ell_{A}^{e_{A}} \ell_{B}^{e_{B}} \ell_{C}^{e_{C}} f \pm 1\)로 서로 다른 세개의 소수를 사용하므로 5차 아이소제니가 필요하며, BSIDH [5]를 구성할 때에는 작은 차수의 아이소제니를 적용 가능할수록 효율적이기 때문에, 낮은 홀수 차수 아이소제니 연산에 대한 효율성 분석도 필요하다. 따라서 본 논문에서는 홀수 차수 아이소제니 연산을 사용했을 때, 어느 형태의 타원곡선이 효율적인지 분석한다. 본 논문의 분석결과 몽고메리와 허프 곡선에서 연산량은 거의 유사하며, 몽고메리와 허프 곡선이 에드워드 곡선보다 0.8% 빠르다는 것을 알 수 있다.

본 논문의 구성은 다음과 같다. 먼저 2장에서는 본 논문 구현 시 사용될 프로토콜인 SIDH에 대해 소개하고, 본 논문에서 사용될 타원곡선 형태에 대해 알아본다. 3장에서는 2장에서 소개된 타원곡선에 대해 SIDH를 구성하는 각 함수의 연산량을 분석한다. 4장에서는 홀수 차수 아이소제니만 사용한 SIDH 구현을 위한 파라미터 세팅 및 구현결과를 제시하고, 5장에서 본 논문을 마무리한다.

Ⅱ. 배경 지식

본 장에서는 구현에서 사용될 타원곡선에 대한 소개를 간략하게 한 뒤, SIDH 프로토콜에 관해 설명한다.

2.1 타원곡선의 종류

2.1.1 몽고메리 곡선과 그 연산

유한체 K 위에서의 몽고메리 곡선은 다음과 같이 정의된다.

Ma,b : by2 = x3 +ax2 +x.       (1)

위 식에서 b(a2 -4)≠0을 만족한다. 특히, b = 1일 경우 간단히 Ma로 나타낸다. SIDH에 사용되는 유한체에서는 몽고메리 곡선 Ma,b는 b = 1인 몽고메리 곡선 Ma로 변환이 가능하므로, 본 논문에서는 Ma인 형태를 사용한다.

Velu의 공식을 이용해서 아이소제니 연산을 수행할 때에는 projective coordinate와 projective curve coefficient를 사용하는 것이 더 효율적으로 연산할 수 있다. \(C \in \bar{K}^{X}\) 에 대하여 a = A/C, b = B/C라 하면, 식 (1)은 다음과 같이 표현할 수 있다.

MA:B:C : By2 = Cx3 +Ax2 +Cx.

Projective coordinate의 경우 다양한 방법이 제안되어있으나, 몽고메리 곡선은 x좌표만 이용해 효율적인 타원곡선 연산이 가능하다는 특징이 있으며, 이에 대한 projective coordinate를 XZ -coordinate라 정의한다 [14]. 다시 말해, 몽고메리 곡선 Ma위의 점 P = (x1,y1)에 대해 함수 x(∙) 를 점에 대한 x좌표를 의미한다. 위의 경우 x(P) = x1을 의미하며, P를 XZ-coordinate로 변환하면 P = (x1: Z1)으로 표현할 수 있다. 이때, x1= x1/Z1을 만족한다.

2.1.2 에드워드 곡선과 그 연산

유한체 K 위에서의 에드워드 곡선은 다음과 같이 정의된다.

Ed : x2 +y2 = 1+dx2y2.       (2)

위 식에서 d≠0,1 이다. 에드워드 곡선에서의 항등원은 (0,1)이며 점 (0,-1)은 위수가 2인 점이다. 또한, (1,0), (-1,0)은 위수가 4인 점이다. Ed위의두 점 P = (x1,y1), Q = (x2, y2)에 대해 P+Q는 다음과 같이 연산한다.

\(P+Q=\left(\frac{x_{1} y_{2}+y_{1} x_{2}}{1+d x_{1} x_{2} y_{1} y_{2}}, \frac{y_{1} y_{2}-x_{1} x_{2}}{1-d x_{1} x_{2} y_{1} y_{2}}\right)\)

타원곡선 점 P에 대한 doubling 연산 2P도 같은 공식을 이용해 구할 수 있다. 마찬가지로 에드워 드 곡선도 projective curve coefficient를 이용해 표현한다. \(C \in \bar{K}^{X}\)에 대하여 d=D/C라 하면, 식 (2)는 다음과 같이 표현할 수 있다.

EC:D : Cx2 + Cy2  = C + Dx2y2

에드워드 곡선에서의 projective coordinate로는 [7]에서 제안된 w-coordinate를 사용한다. [7]에서는 에드워드 위의 점 P = (x,y)에 대하여 w(P) = x2y2로 정의한다. 따라서 P를 WZ -coordinate로 변환하면 P = (W: Z)로 표현할 수 있다.

2.1.3 허프 곡선과 그 연산

허프 곡선은 Joye, Tibouchi, Vergnaud에 의해서 [9]에서 처음으로 제안되었다. 유한체 K위에서의 허프 곡선은 다음과 같이 정의된다.

Ha,b : ax(y2 -1) = by(x2 -1).

위 식에서 a2 ≠b2이고 a,b≠0이다. 허프 곡선에서 항등원은 (0,0)이다. 위의 허프 곡선 Ha,b는 다음과 같이 간단히 나타낼 수 있다.

Hc : cx(y2 -1) = y(x2 -1).       (3)

이 때, c = a/b이고 c≠1을 만족한다. Ha,b위의두 점 P = (x1 ,y1 ), Q = (x2 , y2 )에 대해 P+Q는 다음과 같이 연산한다.

\(\begin{aligned} &P+Q \\ &=\left(\frac{\left(x_{1}+x_{2}\right)\left(1+y_{1} y_{2}\right)}{\left(1+x_{1} x_{2}\right)\left(1-y_{1} y_{2}\right)}, \frac{\left(y_{1}+y_{2}\right)\left(1+x_{1} x_{2}\right)}{\left(1-x_{1} x_{2}\right)\left(1+y_{1} y_{2}\right)}\right) . \end{aligned}\)

에드워드 곡선과 마찬가지로 타원곡선 점 P에 대한 doubling 연산 2P도 같은 공식을 이용해 구할수 있다.

허프 곡선도 projective curve coefficient를 이용해 표현한다. \(C \in \bar{K}^{X}\)에 대하여 C = C/D라 하면, 식 (3)은 다음과 같이 표현할 수 있다.

HC : D : Cx(y2 -1) = Dy(x2 -1).

허프 곡선에서의 projective coordinate로는[8]에서 제안된 w-coordinate를 사용한다. [8]에서는 허프 곡선 위의 점 P = (x,y)에 대하여 w(P) = 1/xy로 정의한다. 따라서 P를 WZ -coordinate로 변환하면 P = (W : Z)로 표현할 수 있다.

2.2 SIDH

서로소인 두 수 ℓA , ℓB에 대해 \(p=\ell_{A}^{e_{A}} \ell_{B}^{e_{B}} f \pm 1\) 형태의 소수를 생성한다. 이때, \(\ell_{A}^{e_{A}} \cong \ell_{B}^{e_{B}}\)가 되도록 eA , eB를 선택한다. 이 소수를 사용해 형성한 유한체 Fp2 위에서 위수가 \(\left(\ell_{A}^{e_{A}} \ell_{B}^{e_{B}} f\right)^{2}\) supersingular 타원곡선 E를 선택한다. 이 경우, ℓ∈{ℓA ,ℓB}와 e∈{eA ,eB} 에 대해 ℓe-torsion subgroup이 생성되며, 기저 {PA ,QA}와 {PB ,QB}를 각각 \(\ell_{A}^{e_{A}}\)\(\ell_{B}^{e_{B}}\) torsion subgroup에서 선택한다.

Alice와 Bob이 서로 비밀 공유키를 교환한다고 가정하자. {PA ,QA}를 Alice의 기저라 하고{PB ,QB}를 Bob의 기저라 하자. 키 생성 단계에서, Alice는 동시에 ℓA로 나누어지지 않는 수 mA , nA\(Z / \ell_{A}^{e_{A}}Z\)에서 선택하고 부분군 〈RA〉=〈[mA]PA+[nA]QA〉를 계산한다. 그 뒤 Velu의 공식을 이용해서 〈RA〉를 커널로 하는 아이소제니 ΦA: E→EA를 연산한다. 여기에서 EA = E/〈RA〉를 만족한다. Alice는 이 아이소제니를 이용해서 연산한 (EA, ΦA(PB), ΦB(QB))를 Bob에게 전달한다. Bob도 비슷한 과정을 통해 Alice에게 (EBB(PA), ΦB(QA))를 연산해 전달한다.

Share secret key를 연산하는 과정에서는, Alice는 Bob에게 받은 점을 이용하여 부분그룹 \(\left\langle R_{A}{ }^{\prime}\right\rangle=\left\langle\left[m_{A}\right] \phi_{B}\left(P_{A}\right)+\left[n_{A}\right] \phi_{B}\left(Q_{A}\right)\right\rangle\)를 연산한 뒤, Velu의 공식을 이용하여 \(\left\langle R_{A}{ }^{\prime}\right\rangle\)를 커널로 하는 아이소제니로 생성된 타원곡선 \(E_{A B}=E_{B} /\left\langle R_{A}^{\prime}\right\rangle\)를 연산한다. Bob도 Alice와 비슷한 과정을 통해 \(E_{B A}=E_{A} /\left\langle R_{B}^{\prime}\right\rangle\)를 생성한다. Alice와 Bob 사이의 비밀 공유키는 생성된 타원곡선의 j-invariant 인 j(EAB) = j(EBA)이다.

Ⅲ. Building blocks for SIDH

본 장에서는 SIDH를 구성하는 각 단위연산에 대해 아이소제니 구현 환경에서의 연산량을 알아본다. 다시 말해, projective coordinate와 projective curve coefficient를 사용했을 때의 연산량을 비교하며, 몽고메리의 경우 XZ-coordinate, 에드워드와 허프 곡선의 경우 WZ-coordinate를 사용한다. 또한, 본 논문에서는 홀수인 소수 — 구체적으로는 3, 5 — 만 사용한 구현이 목표이기 때문에, doubling 연산, differential addition, 3차, 5차 아이소제니 연산에 대해 알아본다. 이 장 이후로 M은 유한체 위에서의 곱셈연산, S는 유한체 위에서의 제곱연산 I는 유한체 위에서의 역원연산을 의미한다.

3.1 몽고메리 곡선

3.1.1 타원곡선 연산

• Doubling: 몽고메리 곡선 Ma에서의 점 P = (X : Z)에 대해 doubling[2]P = (X ′ : Z ′)는 다음과 같이 연산된다.

X ′ = C(X+Z)2(X-Z)2 ,

Z ′ = 4XZ(C(X+Z)2+(A-2C)XZ).

Doubling의 연산량은 4M+2S이다.

• Differential addition: 몽고메리 곡선 위의 두 점 P1 = (X1: Z1 ) P2 = (X2: Z2 )에 대해 x0 = x(P1 -P2 ), x3 = x(P1 +P2 )라 하고, x0 = X0/Z0, x3 = X3/Z3라 하면, 다음과 같이 계산할 수 있다.

X3 = Z0(X1X2 -Z1Z2)2 ,

Z3 = X0(X1Z2 -Z1X2)2 .

Doubling과 differential addition을 동시에 수행하는데 연산량은 6M+4S이다.

3.1.2 아이소제니 연산

• 3-isogeny: 몽고메리 곡선 Ma위의 위수가 3인 점을 P = (X3: Z3 )라 하자. 〈P〉를 커널로 하는 3-isogeny Φ를 Φ : Ma→Ma′ = Ma/〈P〉라 정의하면, Ma의 다른 점 Q = (X : Z)에 대해서 x(Φ(Q))= (X ′ : Z ′)는 다음과 같다.

X ′ = X(XX3 -ZZ3)2 ,

Z ′ = Z(XZ3 -ZX3)2 .

이때 타원곡선의 계수는 이후에 진행되는 tripling 연산을 효율적으로 하기 위해 a 대신에 (a-2)/4를 연산하며, 이에 대한 projective coefficient (A24: C24 ) = (A-2C : 4C)는 다음과 같이 연산한다.

A24 = (X3 +Z3)(Z3 -3X3)3 ,

C24 = 16X3Z33.

3-isogeny의 함수값 (X ′ : Z ′)을 계산하는 과정의 연산량은 4M+2S이며, 이미지 타원곡선의 계수 (A24: C24 )를 계산하는 과정은 2M+3S의 연산량이 필요하다.

• ℓ-isogeny: 몽고메리 곡선 Ma위의 위수가 ℓ = 2s +1인 점을 P = (W1: Z1 )라하고[i]P = (Xi : Zi)라 하자. 〈P〉를 커널로 하는 ℓ-isogeny Φ를 Φ : Ma→Ma′ = Ma/〈P〉라 정의하면, Ma의 다른 점 Q = (X : Z)에 대해서 x(Φ(Q))= (X ′ : Z ′)는 다음과 같다.

\(\begin{aligned} X^{\prime} &=X \prod_{i=1}^{s}\left(X X_{i}-Z Z_{i}\right)^{2} \\ Z^{\prime} &=Z \prod_{i=1}^{s}\left(X Z_{i}-Z X_{i}\right)^{2} \end{aligned}\)

이 때 a′는 다음과 같이 정의한다.

\(a^{\prime}=(6 \tilde{\sigma}-6 \sigma+a) \pi^{2}\)       (4)

위 식에서

\(\sigma=\sum_{i=1}^{s} x\left(P_{i}\right), \tilde{\sigma}=\sum_{i=1}^{s} 1 / x\left(P_{i}\right), \pi=\prod_{i=1}^{s} x\left(P_{i}\right)\)

로 정의한다.

ℓ-isogeny의 함수값 (W ′ : Z ′)을 계산하는 과정의 연산량은 4sM+2S이며, 이미지 타원곡선의 계수 a′에 대해서 a′ = A′/C ′라 할 때, (A ′ : C ′)를 계산하는 과정은 (6s-2)M+3S의 연산량이 필요하다.

3.2 에드워드 곡선

3.2.1 타원곡선 연산

• Doubling: 에드워드 곡선 Ed에서의 점 P = (W: Z)에 대해 doubling [2]P = (W ′ : Z ′)는 다음과 같이 연산된다.

W ′ = 4WZ(D(W+Z)2 -4CWZ),

Z ′ = D(W+Z)2(W-Z)2 .

Doubling의 연산량은 4M+2S이다

• Differential addition: 에드워드 곡선 위의 두 점 P1 = (W1: Z1 ) P2 = (W2: Z2 )에 대해 w0 = w(P1 -P2 ), w3 = w(P1 +P2 )라 하고, w0 = W0/Z0, w3 = W3/Z3라 하면, 다음과 같이 계산할 수 있다.

W3 = Z0(W1Z2 -W2Z1)2 ,

Z3 = W0(W1W2 -Z1Z2)2 .

Doubling과 differential addition을 동시에 수행하는데 연산량은 6M+4S이다.

3.2.2 아이소제니 연산

• 3-isogeny: 에드워드 곡선 Ed위의 위수가 3인 점을 P = (W3: Z3 )라 하자. 〈P〉를 커널로 하는 3-isogeny Φ를 Φ : Ed→Ed′ = Ed/〈P〉라 정의하면, Ed의 다른 점 Q = (W: Z)에 대해서 w(Φ(Q))= (W ′ : Z ′)는 다음과 같다.

W ′ = W(WZ3 -ZW3)2 ,

Z ′ = Z(WW3 -ZZ3)2 .

이때, d′ = D ′/C ′는 다음과 같이 정의한다.

\(\begin{aligned} &D^{\prime}=-16 W_{3}^{2} Z_{3} \\ &C^{\prime}=\left(W_{3}^{2}-2 W_{3} Z_{3}-3 Z_{3}^{2}\right)\left(W_{3}^{2}-6 W_{3} Z_{3}+9 Z_{3}\right)^{2} \end{aligned}\)

3-isogeny의 함수값 (W ′ : Z ′)을 계산하는 과정의 연산량은 4M+2S이며, 이미지 타원곡선의 계수 (D ′ : C ′)를 계산하는 과정은 2M+3S의 연산량이 필요하다.

• ℓ-isogeny: 에드워드 곡선 Ed위의 위수가 ℓ = 2s +1인 점을 P = (W1: Z1 )라 하고[i]P = (Wi : Zi)라 하자. 〈P〉를 커널로 하는 ℓ-isogeny Φ를 Φ : Ed→Ed′ = Ed/〈P〉라 정의하면, Ed의 다른 점 Q = (W: Z)에 대해서 w(Φ(Q))= (W ′ : Z ′)는 다음과 같다.

\(\begin{aligned} W^{\prime} &=W \prod_{i=1}^{s}\left(W Z_{i}-Z W_{i}\right)^{2} \\ Z^{\prime} &=Z \prod_{i=1}^{s}\left(W W_{i}-Z Z_{i}\right)^{2} \end{aligned}\)

이때 d′ = D ′/C ′는 다음과 같이 정의한다.

\(\begin{gathered} D^{\prime}=D^{\ell} \prod_{i=1}^{s}\left(W_{i}+Z_{i}\right)^{8} \\ C^{\prime}=C^{\ell} \prod_{i=1}^{s}\left(2 Z_{i}\right)^{8} \end{gathered}\)

ℓ-isogeny의 함수값 (W ′ : Z ′)을 계산하는 과정의 연산량은 4sM+2S이며, 이미지 타원곡선의 계수(C ′ : D ′)를 계산하는 과정은 2sM+6S+2w(ℓ)의 연산량이 필요하다. 여기에서 w(ℓ) = (h-1)M+(t-1)S를 의미하며, h는 ℓ의 hamming weight를 의미하고 t는 ℓ의 비트길이를 의미한다.

3.3 허프 곡선

3.3.1 타원곡선 연산

• Doubling: 허프 곡선 Hc에서의 점 P = (W: Z)에 대해 doubling[2]P = (W ′ : Z ′)는 다음과 같이 연산된다.

\(\begin{gathered} W^{\prime}=\hat{D}(W+Z)^{2}(W-Z)^{2} \\ Z^{\prime}=4 W Z\left(\hat{D}(W+Z)^{2}+\widehat{C} \cdot 4 W Z\right) \end{gathered}\)

위 식에서 \(\hat{c}=\hat{C} / \hat{D}\)이고, \(\hat{c}=\frac{1}{4}\left(c+\frac{1}{c}-2\right)\)를 만족한다. \(\hat{C}\), \(\hat{D}\)가 주어졌을 때, Doubling의 연산량은 4M+2S이다.

• Differential addition: 허프 곡선 위의 두 점 P1 = (W1: Z1 ) P2 = (W2: Z2 )에 대해 w0 = w(P1 -P2 ), w3 = w(P1 +P2 )라 하고, w0 = W0/Z0, w3 = W3/Z3라 하면, 다음과 같이 계산할 수 있다.

W3 = Z0 (W1W2 -Z1Z2)2 ,

Z3 = W0(W1Z2 -Z1W2)2 .

Doubling과 differential addition을 동시에 수행하는데 연산량은 6M+4S이다.

3.3.2 아이소제니 연산

• 3-isogeny: 허프 곡선 Hc위의 위수가 3인 점을 P = (W3: Z3 )라 하자. 〈P〉를 커널로 하는 3-isogeny Φ를 Φ : Hc→Hc′ = Hc/〈P〉라 정의하면, Hc의 다른 점 Q = (W: Z)에 대해서 w(Φ(Q))= (W ′ : Z ′)는 다음과 같다.

W ′ = W(WW3 -ZZ3)2 ,

Z ′ = Z(WZ3 -ZW3)2 .

이때, c′ = C ′/D ′는 다음과 같이 정의한다.

C ′ = C(DZ3 +CW3)2 ,

D ′ = D(CZ3 +DW3)2 .

3-isogeny 연산 후 얻은 타원곡선에서 tripling 연산을 쉽게 하도록 c′대신에 \(\hat{c}=\frac{1}{4}\left(c^{\prime}+\frac{1}{c^{\prime}}-2\right)=\hat{C} / \hat{D}\)를 연산하는 것이 더 효율적이며, 이떄 \(\hat{C}\), \(\hat{D}\)는 다음과 같이 연산한다.

\(\begin{gathered} \hat{C}=\left(W_{3}+Z_{3}\right)\left(Z_{3}-3 W_{3}\right)^{3} \\ \hat{D}=16 W_{3} Z_{3}^{3} \end{gathered}\)

3-isogeny의 함수값 (W ′ : Z ′)을 계산하는 과정의 연산량은 4M+2S이며, 이미지 타원곡선의 계수 (\(\hat{C}\):\(\hat{D}\))를 계산하는 과정은 2M+3S의 연산량이 필요하다.

• ℓ-isogeny: 허프 곡선 Hc위의 위수가 ℓ = 2s +1인 점을 P = (W1: Z1 )라하고[i]P = (Wi : Zi)라 하자. 〈P〉를 커널로 하는 ℓ-isogeny Φ를 Φ : Hc→Hc′ = Hc/〈P〉라 정의하면, Hc의 다른 점 Q = (W: Z)에 대해서 w(Φ(Q))= (W ′ : Z ′)는 다음과 같다.

\(\begin{aligned} W^{\prime} &=W \prod_{i=1}^{s}\left(W W_{i}-Z Z_{i}\right)^{2} \\ Z^{\prime} &=Z \prod_{i=1}^{s}\left(W Z_{i}-Z W_{i}\right)^{2} \end{aligned}\)

이때, c′ = C ′/D ′는 다음과 같이 정의한다.

\(\begin{aligned} &C^{\prime}=C \prod_{i=1}^{s}\left(D Z_{i}+C W_{i}\right)^{2} \\ &D^{\prime}=D \prod_{i=1}^{s}\left(C Z_{i}+D W_{i}\right)^{2} . \end{aligned}\)

ℓ-isogeny의 함수값 (W ′ : Z ′)을 계산하는 과정의 연산량은 4sM+2S이며, 이미지 타원곡선의 계수(C ′ : D ′)를 계산하는 과정은 4sM+2S의 연산량이 필요하다.

Remark 1. ℓ-isogeny도 마찬가지로 이미지 타원곡선 위에서 효율적인 타원곡선 연산을 수행하기 위해서는 c′가 아닌 \(\hat{c}=\frac{1}{4}\left(c^{\prime}+\frac{1}{c^{\prime}}-2\right)\)를 구하는 것이 더 편리하다. 하지만 3-isogeny 공식처럼 타원 곡선의 계수 구하는 공식을 커널의 원소로만 표현하기 위해서는 division polynomial이 필요하며, 5차 이상 아이소제니에 대해서 이를 이용해 식을 일반화하기란 어렵다. 따라서 일반적으로 5차 이상 아이소제니 연산을 수행할 때에는 c′를 계산한 뒤에, \(\hat{C}\)로 변환하는 과정을 거쳐야 한다. 다시 말해 c′ = C ′/D ′에 대해 \(\hat{C}=\left(C^{\prime}-D^{\prime}\right)^{2}\), \(\hat{D}=4 C^{\prime} D^{\prime}\)를 구하는 추가적인 연산이 필요하며, 이 연산은 2S가 필요하다.

다음 Table 1.은 몽고메리, 에드워드, 허프 곡선 위에서 SIDH를 구성하는 함수들의 연산량을 비교한 표이다. 아래 표에서 Mont.는 몽고메리 곡선을, Edwards는 에드워드 곡선, Huff는 허프 곡선을 의미한다. DBLADD는 doubling 연산과 differential addition 연산을 의미하고 DBL은 doubling 연산을 의미한다. 3-isogeny는 3차 아이소제니 함수값을 계산하는 연산과 타원곡선 계수를 계산하는 연산을 합한 연산량을 의미한다. 또한, ℓ-eval은 ℓ차 아이소제니 함수값을 계산하는 연산량을 의미하며, ℓ-coeff은 ℓ차 아이소제니의 이미지 곡선의 계수를 복원하는 연산량을 의미한다. coeffTrans는 타원곡선의 계수를 타원곡선 연산하기 쉬운 형태로 변환하는 과정을 의미하며, 이는 허프 곡선에서만 사용됨을 알 수 있다. 특히, 본 논문에서는 3차, 5차 아이소제니를 활용해 구현을 할 것이므로, 5차 아이소제니에서의 함수값을 계산하는 연산량과 타원곡선의 계수를 복원하는 연산량을 비교하였다.

Table 1. Computational cost of the building blocks of SIDH

JBBHCB_2021_v31n1_19_t0001.png 이미지

위의 표에서 확인할 수 있듯이, 몽고메리, 에드워드, 허프 곡선의 타원곡선 연산량 및 3차 아이소제니 연산량이 같음을 알 수 있다. 하지만 아이소제니 차수가 올라갈수록 몽고메리보다 에드워드나 허프 곡선이 더 효율적임을 알 수 있다. 다음 장에서는 추가적인 최적화 기법을 사용해서, 3차 5차 아이소제니를 이용해 SIDH를 구현하는 방법과 구현결과를 제시한다.

Ⅳ. 구현결과

본 장에서는 3차, 5차 아이소제니를 활용한 SIDH를 구현결과를 제시하였는데, 해당 차수를 선정한 이유는 다음과 같다.

• 에드워드 곡선의 w-coordinate의 경우 4차 compression 함수로, 효율적인 4차 아이소제니를 구성하기가 힘들다. 허프 곡선의 경우 3차, 4차 아이소제니를 사용한 SIDH 구현은 이미 [8]에 제시되어있다.

• 따라서 현재 아이소제니 기반 암호 구현에 언급되고 있는 몽고메리, 에드워드, 허프 곡선을 동일선상에서 비교하기 위해서는 홀수 차수 아이소제니를 활용하는 것이 좋다.

• 특히 허프 곡선의 경우 5차 아이소제니 연산 시타원곡선 계수 변환이 중요한 점으로 작용할 수 있으며, 몽고메리 곡선의 경우 SIDH 구현에서는 효율적인 타원곡선 계수 변환 방법이 존재하기 때문에, 이를 적용해 구현 성능을 알아보는데 5차 아이소제니가 적합하다고 볼 수 있다.

따라서 본 장에서는 홀수 차수 아이소제니를 이용한 SIDH를 구현하기 위해 설정한 파라미터와 추가적인 구현 기법, 그리고 구현결과를 제시한다.

4.1 파라미터 설정

홀수 차수 아이소제니를 이용한 SIDH 구현을 위해 사용한 621-비트 소수는 다음과 같다.

p = 26731755119 -1.

위의 소수에 대해 3175 ≅2277.368이며 5119 ≅ 2276.309로, SIDH 프로토콜이 요구하는 \(\ell_{A}^{e_{A}} \cong \ell_{B}^{e_{B}}\)

를 만족시킴을 알 수 있다. 해당 소수를 사용할 경우 고전 컴퓨팅 환경에서 138비트의 보안강도를 가진다. 위의 소수를 이용한 유한체 Fp2 = Fp (i) 위에 정의된 다음 몽고메리, 에드워드, 허프 곡선을 이용하였다.

M : y2 = x3 +x,

E-1: x2 +y = 1-x2y2 ,

Hi : ix(y2 -1) = y(x2 -1).

Fp2에서 SIDH 구현을 위해서는 generator point PA ,QA, PB ,QB의 선택이 필요하다. 이때, PA ,QA∈E[3175]에 있으며, 두 점 모두 위수가 3175이다. 또한, PB ,QB∈E[5119]이며 두 점 모두 위수가 5119이다. 해당 점을 선택하는 방법은 SIKE specification [1]에 나와 있는 방법을 따랐으며, 몽고메리 곡선에서 generator point를 선택하고, 해당 점을 에드워드 곡선으로 이동해서 w -coordinate로 변환해서 사용하였다. 허프 곡선의 경우 다음 식과 같은 동형인 Weierstrass curve를 사용하였다.

W : y2 = x3 -x.

위의 곡선에서 generator point를 설정한 이후 허프 곡선 위의 점으로 변환하여서 사용하였다.

4.2 몽고메리 타원곡선 계수 복원 기법

에드워드와 허프 곡선의 경우 타원곡선 계수를 계산하는 방법이 비교적 최적화 되어있으나, 몽고메리 곡선의 경우 5차 이상의 아이소제니의 경우 타원곡선 계수를 복원하는 방법이 다소 복잡하다. 본 논문에서는 5-isogeny 구현 시, [4]에서 제시된 2-torsion을 이용하는 방법으로 타원곡선의 계수를 변경하였고 이 방법에 대해 상세히 설명한다.

몽고메리 곡선 Ma는 항상 다음과 같은 형태의 위수가 2인 점 (2-torsion point)를 가진다.

P0 = (0,0),Pα (α,0),P1/α (1/α,0).

이를 projective XZ-coordinate로 변환하면 다음과 같이 표현된다.

P0 = (0 : 1),Pα = (Xa : Za),P1/α = (Za : Xa).

또한, Pα는 Ma위의 점이므로, a =-(α2 +1)/α를 만족한다. 몽고메리 곡선에서는 연산의 편의를 위해 a대신 (a-2)/4를 사용하는데, 이를 다시 projective coordinate로 표현하면

(A-2C : 4C) = ((Xα +Zα )2:-4XαZα ).

로 나타낼 수 있으며, 해당 연산은 2S가 소요된다. 다시 말해 2-torsion point P에 대해서 Φ(P)가 주어지면, 이를 이용하여 이미지 곡선의 계수를 변경할 수 있음을 알 수 있다. 특히, 본 논문에서는 홀수 차수 아이소제니를 사용하기 때문에 2-torsion point P에 대한 아이소제니 연산을 여러 번 하여도, P의 위수는 변하지 않고 2로 유지된다. 따라서, SIDH 프로토콜을 처음 수행할 때 base curve에서의 2-torsion point를 공개값으로 가지고 있고, 이에 대한 아이소제니 연산과 계수 복원연산을 수행하면 이미지 곡선의 계수를 구할 수 있다.

그 후 상대방이 연산한 타원곡선에서 공유할 비밀 값을 연산하는 과정에서는 몽고메리 곡선에서 이차방정식을 푸는 연산을 이용해 2-torsion point를 직접 계산한 다음, 해당 점을 이용해서 다시 타원곡선 계수 복원을 진행한다. 이차방정식을 풀 때는 square-root 는 한 번 진행되고, p≡3mod4를 만족하므로 간단한 지수연산으로 수행할 수 있다. 또한, 해당 연산은 Bob이 shared secret을 구할 때만 일어나므로, 전체 SIDH 과정에서 한 번만 수행된다.

다음 표는, 몽고메리 곡선에서 식 (4)를 사용해서 5차 아이소제니의 이미지 곡선의 계수를 복원하는 경우와 2-torsion point를 이용했을 때의 연산량을 비교한다. Table 2.에서 Original은 식 (4)를 이용한 방법을, 2-torsion은 2-torsion을 이용해서 복원한 방법을 의미한다.

Table 2. Comparison of the computational costs for recovering the coefficient of the image curve

JBBHCB_2021_v31n1_19_t0003.png 이미지

Remark 2. SIDH 기반 암호의 경우, 상대방의 generator Pi,Qii∈{A,B}와 자신의 개인키 n을 이용해 Pi +nQi를 연산하는 부분이 있다. 해당 부분은 몽고메리 래더 (Montgomery ladder) 방식을 이용해 구현하면 효율적으로 구현할 수 있으며, 따라서 SIDH는 Pi ,Qi 뿐만 아니라 Ri = Pi -Qi도 generator point로 여겨 연산한다. [4]에서는 두 점 P, Q와 그의 차 P-Q를 이용하면 몽고메리 곡선 계수 a를 복원하는 방법에 대해서 제안하며, 해당 방법을 differential method로 명시하겠다.

Differential method를 이용해 계수를 복원하는 방법은 8M+5S의 연산량이 소요된다. 2-torsion point와 비교해서 squaring이 한 번 밖에 더 들지 않으므로, 이 방법이 2-torsion을 이용하는 방법보다 효율적이라고 생각할 수 있다. 하지만 타원곡선의 계수를 복원하는 과정은 매 아이소제니 연산시 일어나야 하므로, 논문에서 제시한 파라미터상으로는 2-torsion 방법보다 119 × 2 = 238번의 Fp2에서의 제곱연산이 더 필요하다는 것을 의미한다. (public key generation에서 119번, shared secret key generation에서 119번). 따라서 5-isogeny 연산 시에는 2-torsion point 방법이 더 효율적임을 알 수 있다. 하지만 5차 이상의 아이소제니에서는 differential method가 더 효율적이라는 것을 알 수 있다.

4.3 실험 결과

본 장에서는 홀수 차수를 이용한 SIDH를 몽고메리, 에드워드, 허크 곡선에서 구현한 결과를 제시한다. 측정에 사용한 CPU는 3.40GHz의 동작 주파수를 가지는 Intel Core i7-6700를 사용했으며, Ubuntu 18.04.2 LTS 운영체제상에서 clang-6.0.0 컴파일러를 이용했다. Table 3.에서 A Keygen은 Alice의 공개키 생성을 의미하며, B Keygen은 Bob의 공개키 생성을 의미한다. 마찬가지로 A Keyshare는 Alice의 shared secret key 연산, B Keyshare는 Bob의 shared secret key 연산을 의미한다. 본 실험에서 Alice가 3차 아이소제니를 사용하고 Bob이 5차 아이소제니를 사용하는 것으로 한다.

Table 3. Performance result of SIDH implementation. The results were rounded to the nearest 103 clock cycles

JBBHCB_2021_v31n1_19_t0002.png 이미지

Table 3.에서 확인할 수 있듯이, 몽고메리와 허프 곡선에서의 연산량은 유사하며, 에드워드 곡선에서는 몽고메리나 허프 곡선과 비교했을 때 0.8% 느리다는 것을 알 수 있다.

먼저 몽고메리와 허프 곡선을 결과를 분석해본다. 3차 아이소제니를 사용하는 Alice 쪽은 몽고메리와 허프 곡선에서 SIDH를 구성하는 모든 단위연산이 동일하여 연산시간이 같다는 것을 쉽게 알 수 있다. 몽고메리와 허프 곡선의 연산량이 달라지는 부분은 Bob 쪽에서 5차 아이소제니의 이미지 곡선의 계수를 복원하는 부분이다.

허프 곡선의 경우 5차 아이소제니의 이미지 곡선 계수를 복원하는데 소요되는 연산량은 8M+2S이다. 하지만 생성된 계수를 연산에 효율적으로 바꾸기 위해서는 Table 1.의 coeffTrans 연산이 필요하며, 따라서 총 8M+4S의 연산량이 필요하다. 몽고메리 곡선의 경우 이미지 곡선의 계수를 복원하는데 2-torsion 방법이 사용될 수 있으며, 이 방법은 8M+4S의 연산량이 필요하다. 따라서, B Keygen에서도 몽고메리 곡선과 허프 곡선 사이에서 SIDH 구성하는 모든 단위 연산량이 일치하여, 사이클 수도 비슷하게 나온다는 사실을 확인할 수 있다.

B Keyshare에서 몽고메리와 허프 곡선을 비교해보면, 몽고메리의 경우 Alice로부터 받은 타원곡선의 계수를 이용하여 2-torsion 점을 얻어내는 과정이 필요하다. 주어진 타원곡선의 계수를 a라 하면, 이는 x2 +ax +1의 근을 구하는 과정이고, Fp2에서 연산이 진행되기 때문에 해는 존재한다. 해당 해를 구하는데 크게 1번의 제곱근 연산이 필요하다. 한편, 허프 곡선의 경우 Hc에 대해서 \(\hat{c}=c+1 / c+2\) 가 전달된다. Alice의 경우 \(\hat{c}\)를 이용해서 바로 타원곡선 및 아이소제니 연산에 사용할 수 있지만, Bob의 경우 이미지 곡선의 계수를 복원하기 위해서는 c자체가 필요하다. c도 이차방정식의 해를 구하는 과정이 필요하고, 크게 1번의 제곱근 연산이 필요하다. 결론적으로 3차, 5차 아이소제니 사용에 있어서 SIDH 를 구성하는 연산량이 몽고메리와 허프 모두 동일하기 때문에, 두 구현결과는 비슷한 사이클 수를 가진다는 것을 해석할 수 있다.

한편, Table 3.에서 살펴보면 몽고메리와 허프 곡선을 이용해 SIDH를 구현한 것보다 에드워드 곡선에서 0.8% 느리다는 것을 확인할 수 있다. Alice 쪽에서의 3차 아이소제니 연산의 경우 몽고메리, 에드워드, 허프 곡선에서의 공식이 동일하여 연산량의 차이가 거의 발생하지 않는다. 하지만 Bob 쪽에서 5차 아이소제니 연산에서는 이미지 곡선의 계수를 복원하는 과정에서 연산량의 차이가 발생한다. 에드워드 곡선에서는 이미지 타원곡선의 계수를 복원하는데 6M+10S의 연산량이 필요하며, 1S≅0.8M 인 부분을 생각하면, 해당 연산량은 14M이 필요하다. 반면 허프 곡선이나 몽고메리 곡선에서는 8M+4S가 필요하며, 이는 약 11.2M에 해당한다. Bob에서 이미지 타원곡선의 계수를 복원하는 부분이 총 119 × 2 = 238번 (public key generation에서 119번, shared secret key generation에서 119번)이 필요하기 때문에, 이는 에드워드 곡선에서의 속도 저하를 가져온다.

V. 결론

본 논문에서는 SIDH 프로토콜을 홀수 차수 아이소제니만 이용하여 구현할 경우 몽고메리, 에드워드, 허프 곡선 중 어느 곡선이 효율적인지에 관해 실험을 진행하였다. 이와 관련하여, 먼저 각 곡선에서 SIDH 구현을 위한 단위연산의 연산량을 비교하였으며 홀수 차수 사용을 위한 유한체를 새롭게 정의하였다. 본 논문의 결과 에드워드 곡선이 몽고메리나 허프 곡선보다 0.8% 느리다는 것을 알 수 있으며, 몽고메리와 허프 곡선에서는 속도가 유사하다는 점을 알 수 있다. 특히, differential method를 사용할 경우 몽고메리와 허프 곡선에서 타원곡선 계수 복원은 아이소제니의 차수와 관계없이 일정하므로, SIDH 기반 암호에서는 차수가 올라갈수록 몽고메리, 허프 곡선과 에드워드 곡선의 성능 차이가 발생할 수 있다.

* 이 논문은 정부(과학기술정보통신부)의 재원으로 한국연구재단의 지원을 받아 수행된 연구임(No. NRF-2017R1A2 B4011599)

References

  1. R. Azarderakhsh et al. "Supersingular isogeny key encapsulation", submission to the NIST post-quantum standardization project, 2017
  2. R. Azarderakhsh et al. "Practical supersingular isogeny group key aggrement," IACR Cryptology ePrint Archive, 2019:330, 2019
  3. J. Couveignes, "Hard homogeneous spaces," IACR Cryptology ePrint Archive, 2006:291, 2006
  4. C. Costello and H. Hisil, "A simple and compact algorithm for SIDH with arbitrary degree isogenies," Advances in Cryptology, ASIACRYPT'17, LNCS 10625, pp. 303-329, 2017
  5. Craig Costello "B-SIDH supersingular isogeny Diffie-Hellman using twisted torsion," Advances in Cryptology, ASIACRYPT'20, LNCS 12492, pp. 440-463,2020
  6. R. Drylo et al. "Efficient Montgomerylike formulas for general Huff's and Huff's elliptic curves and their applications to the isogeny-based cryptography," IACR Cryptology ePrint Archive, 2020:526, 2020
  7. R. Farashahi et al. "Differential addition on twisted Edwards curves," ACISP'17, LNCS 10343, pp. 366-378, 2017
  8. Y. Huang et al, "Optimized arithmetic operations for isogeny-based cryptography on Huff curves," ACISP'20, LNCS 12248, pp. 23-40, 2020
  9. M. Joye et al,"Huff's model for elliptic curves," International Algorithmic Number Theory Symposium, ANTS'10, pp. 234-250, 2010
  10. D. Jao, L. De Feo "Towards quantum-resistant cryptosystems from supersingular elliptic curve isogenies," PQCrypto'11, LNCS 7071, pp. 19-34, 2011
  11. S. Kim et al. "Optimized method for computing odd-degree isogenies on Edwards curves," Advances in Cryptology, ASIACRYPT'19, LNCS 11922, pp. 273-292, 2019
  12. S. Kim et al. "New hybrid method for isogeny-based cryptosystems using Edwards curves," IEEE transactions on Information Theory, vol. 66, no. 3, pp. 1934-1943, 2020 https://doi.org/10.1109/tit.2019.2938984
  13. M. Meyer et al. "On hybrid SIDH schemes using Edwards and Montgomery curve arithmetic," IACR Cryptology ePrint Archive, 2017:1213, 2017
  14. P. Montgomery, 'Speeding the pollard and elliptic curve methods of factorization," Mathematics of computation, vol. 48, no. 177, pp. 243-264, 1971 https://doi.org/10.1090/S0025-5718-1987-0866113-7
  15. D. Moody and D. Shumow, "Analogues of Velu's formula for isogenies on alternate models of elliptic curves," Mathematics of Computations, vol. 85, no. 300, pp. 1929-1951, 2016
  16. P. W. Shor, "Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer," SIAM review, vol. 41 no. 2, pp. 303-332, 1999. https://doi.org/10.1137/S0036144598347011