https://keymoon.hatenablog.com/entry/2024/12/05/032548

この記事はCTF Advent Calender 2024の4日目の記事です。

昨日の記事はn01e0さんのSECCON 13 Quals - Jumpでした。作問ミスって怖いですね。カツカツスケジュールになってしまうのも作問ミスも他人事ではないので、襟を正してレビュー等をしようと思います。

この記事は元々はCryptoワザップと題して「400-bit程度の合成数までならノートPCで素因数分解できますよ〜、512-bit程度の合成数ならば5000円と数時間くらいで素因数分解できますよ〜」みたいなことを書こうとしていた記事でした。しかし、その一環で楕円曲線法について調べ出したらRabbit Holeに落ちてしまったので、楕円曲線法に限定した記事として公開することにしました。

本記事では、楕円曲線法のアルゴリズムを解説した後、楕円曲線法に適用される高速化のテクニックを紹介します。実装やソフトウェアの使用方法などには触れません。

記号・用語:

$n$ : 素因数分解する数
$p$ : 求めたい素因数
$E(K)$ : $K$ 上の楕円曲線
$P_x$ : 楕円曲線上の点 $P$ のアフィン座標におけるx（ $P_y$ も同様）
$P_Z$ : 楕円曲線上の点 $P$ の射影座標（後述）におけるZ
$0_E$ : 楕円曲線 $E$ の無限遠点

楕円曲線法の概要

まずは、簡単に楕円曲線法のアルゴリズムについて紹介します。楕円曲線法はp-1法やp+1法の原理を応用したアルゴリズムです。これらのアルゴリズムは $p$ が計算量に大きな影響を与えるという特性があります。これは、計算量が $n$ に大きく影響される二次ふるい法や数体ふるい法とは対照的です。そして、楕円曲線法はこれらのアルゴリズムのなかでも最も高い性能を持つアルゴリズムです。そのため、楕円曲線法は $n$ の大きさに依存する素因数分解アルゴリズムである数体ふるい法などを適用する前処理として用いられます。他にも、数体ふるい法を適用できないような非常に大きな数を素因数分解するのに用いられたりもします。2024年現在、楕円曲線法で見つかった最大の素因数は7^337+1の素因数である274-bitの数です(ECMNET)。

p-1法やp+1法、楕円曲線法のアルゴリズムについては、以下の記事に分かりやすくまとまっています。

wacchoz.hatenablog.com

本記事でも、簡単にですが楕円曲線法自体について説明をします。

楕円曲線法は、 $n$ の非自明な約数をランダムに発見する乱択アルゴリズムです。発見できる約数の大きさや約数を発見できる確率は、アルゴリズムに与えるパラメータによって変動します。以降では、説明の簡略化のためにアルゴリズムが成功した場合は素因数 $p$ を発見するものとします*1。

以下は、高速化が行われていない教科書的な楕円曲線法の大まかな実装です。

def get_M(n, bound):
  m = 1
  # equivalent to: for i in range(1, bound): m = LCM(m, i)
  for p in primes(bound+1): # primes in [0, bound]
    k = floor(log(bound) / log(p)) # k such that p^{k} <= bound < p^{k+1}
    m *= p**k
  return m

def ecm(n, bound):
  E, P0 = get_random_elliptic_curve_and_point(N)
  m = get_M(n, bound)
  try:
    Q = m * P0
  except ModInvException as e:
    return gcd(e.value, n)
  return None

楕円曲線法では、まず $\mathbb{Z}/n\mathbb{Z}$ 上の楕円曲線（のようなもの*2） $E(\mathbb{Z}/n\mathbb{Z})$ をランダムに生成します。アルゴリズムが素因数 $p$ の発見に成功できるとして、 $E(\mathbb{Z}/n\mathbb{Z})$ と同じパラメータの $\mathbb{Z}/p\mathbb{Z}$ 上の楕円曲線 $E(\mathbb{Z}/p\mathbb{Z})$ を考えてみます*3。ここで、 $P+Q=0_{E(\mathbb{Z}/p\mathbb{Z})}$ 、つまり $P=-Q$ になるような点 $P,Q$ を見つけられたとします。このとき、元の $E(\mathbb{Z}/n\mathbb{Z})$ で $P+Q$ はどうなるでしょうか？楕円曲線の性質より、 $P=(P_x, P_y)\equiv (Q_x, -Q_y)=-Q \bmod p$ が成立します。よって、 $P_x\equiv Q_x \bmod p$ です。一方で、楕円曲線の加法公式には $Q_x\neq P_x$ の場合に $L:=\frac{Q_y-P_y}{Q_x-P_x}$ が係数として登場します。ほとんどの場合 $\gcd(Q_x-P_x, n)=p$ なので、この $L$ は計算することができません。よって、 $P+Q$ を満たす値は $E$ 上に存在しないことが分かりました。同時に、 $P+Q$ を計算する過程で素因数 $p$ を発見できていることも分かるかと思います。

ここまでの議論で、 $E(\mathbb{Z}/p\mathbb{Z})$ で $P+Q=0_E(\mathbb{Z}/p\mathbb{Z})$ となるような $P$ と $Q$ を発見できれば、素因数 $p$ を発見することができると分かりました。これは、加算ではなくスカラー倍でも同様の結論になります*4。つまり、 $mP\equiv 0_{E(\mathbb{Z}/p\mathbb{Z})}$ を満たすような正整数 $m$ と $P$ を発見することでも、素因数 $p$ を発見することができるということです。 $mP=0_{E(\mathbb{Z}/p\mathbb{Z})}$ は、 $m$ が $E(\mathbb{Z}/p\mathbb{Z})$ の位数 $\#E(\mathbb{Z}/p\mathbb{Z})$ の倍数であるときかつそのときに限り成立します。ここで、うまいこと $m$ をとれば、 $\#E(\mathbb{Z}/p\mathbb{Z})$ の一定割合は $m$ を割り切ってくれるのではないかというのが楕円曲線法のアイデアです。

$m$ の定め方は、あるパラメータ $B$ を用いて、 $m:=LCM(1, 2, 3, ..., B)$ とします。こうすると、 $\#E(\mathbb{Z}/p\mathbb{Z})$ が $B$ -power-smoothな場合、つまり各素因数毎の積 $p^e$ が全て $B$ 以下の場合にアルゴリズムが成功すると言えます。

$\#E(\mathbb{Z}/p\mathbb{Z})$ は $[p+1-2\sqrt{p}, p+1+2\sqrt{p}]$ をとります*5(Hasse's theorem)。このことを用いて、求めたい $p$ の大きさ毎にパラメータ $B$ をうまく定めて素因数分解を行います。

計算量の解析

$\rm{M}(\log N)$ を $\log_2 N$ -bitの整数の乗算にかかる計算量とします。すると、楕円曲線法の具体的な計算量は $O(\exp( (\sqrt{2}+o(1)) \cdot \sqrt{\log p\log\log p}) \rm{M}(\log N) )$ となります。本項では、これを導出します。この計算量は、先程のecmルーチン一度の実行にかかる計算量ではなく、 $p$ 以下の素因数を見逃す確率が $e^{-1}$ を下回るために必要な計算量であることに注意してください。

まず、先程のecmルーチンに必要な計算量を見積もります。boundを $B$ とおくと、 $m \approx B^{\pi(B)}$ が成立します。そのため、 $mP$ を計算するために必要な計算量は素数定理を用いることで

$\begin{align} O(\log B^{\pi(B)} \rm{M}(\log N)))&=O(\pi(B)\cdot \log B \rm{M}(\log N)))\\ &=O(\frac{B}{\log B}\cdot \log B \rm{M}(\log N)))\\ &=O(B \rm{M}(\log N))) \end{align}$

と分かります。これはecmルーチンで最も支配的な部分であるため、これがルーチンの計算量となります。

次に、ecmルーチンを実行する回数と $B$ の最適な値を見積もります。

まず、関数 $L$ を以下のように定義します。

$L(x) = \exp(\sqrt{\log x \log\log x})$

$x$ 以下の整数が $L(x)^a$ -smooth*6である確率は、 $x\to \infty$ で $L(x)^{-1/(2a)+o(1)}$ であることが知られています(Canfield-Erdős-Pomerance theorem)。これは整数を $[x-2\sqrt{x}, x+2\sqrt{x}]$ から選んだときにも、経験的に同じことが言えます。

$B$ を $L(p)^a$ にして $1/(L(p)^{-1/(2a)+o(1)})$ 回ルーチンを実行することを考えます。これの成功確率は、 $(1-\epsilon)^{1/\epsilon} \to_{\epsilon \to 0} e^{-1}$ であることから $p \to \infty$ で $1-e^{-1}$ に漸近します。なので、このパラメータでと実行回数で要求する成功確率を達成できています。最後に、全体の計算量を考えます。これは、 $1/(L(p)^{-1/(2a)+o(1)})\cdot O(L(p)^a \rm{M}(\log N))=O(L(p)^{1/(2a)+a+o(1)} \rm{M}(\log N))$ となります。この式を最適化する $a$ は $\sqrt{2}$ であるため、これを代入して $O(L(p)^{\sqrt{2}+o(1)} \rm{M}(\log N))$ を得ることができます。

楕円曲線法の高速化

教科書的な実装を再掲します:

def get_m(n, bound):
  m = 1
  # equivalent to: for i in range(1, bound): m = LCM(m, i)
  for p in primes(bound+1): # primes in [0, bound]
    k = floor(log(bound) / log(p)) # k such that p^{k} <= bound < p^{k+1}
    m *= p**k
  return m

def ecm(n, bound):
  E, P0 = get_random_elliptic_curve_and_point(n)
  m = get_m(n, bound)
  try:
    Q = m * P0
  except ModInvException as e:
    return gcd(e.value, n)
  return None

この実装を高速化していきます。

高速化手法1: 射影座標とLucas Chain、基礎算術の高速化

まず、楕円曲線上の点を射影座標で持つことで、楕円曲線上の加減算で除算を行わないようにできます。射影座標を用いることは原著論文(Lenstra 1985)の時点から指摘されていました。射影座標については、以下で詳しく解説されています。

zenn.dev

次に、 $m$ を陽に計算しないようにし、その代わりにループごとに $Q := pQ$ をk度計算することとします。楕円曲線上のスカラー倍は、既知の点の加算か二倍算を繰り返して行われます。射影座標は他の座標の持ち方と比べて加算に多くのステップ数がかかるという欠点がありますが、 $P-Q$ が分かっているならば高速に計算することができます（詳細は20 years of ECMの2節を参照）。そのため、加算は差が既知の点同士のみで行いたいです。このような条件下で、最小ステップ数で $pQ$ を計算するための加法連鎖はLucas Chainとして知られています。Lucas Chainの長さは $O(\log n)$ となるので、加法連鎖の条件に制限をつけることでのペナルティは大きくありません。以下は、 $23Q$ を計算するためのLucas Chainの一例です:

2P := 2*P
3P := 2P + P (差はP)
5P := 2P + 3P (差はP)
7P := 2P + 5P (差は3P)
9P := 2P + 7P (差は5P)
16P := 7P + 9P (差は2P)
23P := 7P + 16P (差は9P)

このようなLucas Chainは、小さい $p$ については前計算することができます。Lucas Chainが既知でない場合も、PRACアルゴリズムというアルゴリズムを用いることで最適な演算回数に近いchainを生成することができます*7。

さらに、これらのアルゴリズム演算の根底にある多倍長演算やmod p上での乗算の高速化も行えます。これには998244353で割ったあまりを求めることを日課にしている勢力にはお馴染みの、モンゴメリ乗算やBarrett Reductionといったテクニックを適用することができます。

natsugiri.hatenablog.com

これらの最適化を実装したコードを以下に示します:

def ecm(n, bound):
  E, P0 = get_random_elliptic_curve_and_point(n)
  Q = P0.to_proj_coords()
  for p in filter(is_prime, range(bound)): # primes in [0, bound)
    k = floor(log(bound) / log(p)) # k such that p^{k} <= bound < p^{k+1}
    p_chain = get_chain(p)
    for i in range(k):
      Q = compute_mul_using_chain(Q, p_chain)
    if gcd(Q.Z, n) != 1: return gcd(Q.Z, n)
  return None

高速化手法2-1: standard continuation

次に、アルゴリズム自体に改良を加えます。新しい処理を今までのアルゴリズムに続く形で実装し、今までのアルゴリズムをstage 1、加えたアルゴリズムをstage 2と呼ぶことにします。また、stage 1で用いていたboundをbound1(あるいはB1)と言い換え、stage 2で用いる新たなboundをbound2(あるいはB2)と呼ぶことにもします。

stage 2では、stage 1で計算された点である $Q:=m\cdot P_0$ に様々な係数を別々で掛け合わせます。具体的には、 $[\rm{B1}, \rm{B2})$ に含まれる素数すべてについて、 $\gcd((pQ)_Z, n)$ を都度確かめます。こうすることで、列挙する素数が足りず「あと一歩」だった $E$ に対しても、素因数分解を行うことができます。これには一見すると嬉しさが殆どないように思えるかもしれませんが、このステップがあることで $log p$ 倍ほどの高速化がもたらされます(Brent 1985)。また、以降に述べるテクニックによってこの計算はより高速に行えるようになります。

def ecm_stage1(n, E, B1): # not optimized
  return get_m(n, B1) * E.random_element()

def ecm_stage2(n, E, Q, B1, B2):
  for p in primes(B1, B2+1): # primes in [B1, B2]
    Q2 = p * Q
    if gcd(Q2.Z, n) != 1: return gcd(Q2.Z, n)
  return None

def ecm(n, B1, B2):
  assert B1 <= B2
  E = get_random_elliptic_curve(n)
  Q = ecm_stage1(n, E, bound1)
  if gcd(Q.Z, n) != 1: return gcd(Q.Z, n)

  return ecm_stage2(n, E, Q, B1, B2)

高速化手法2-2: Montgomery continuation

まず、 $[\rm{B1}, \rm{B2})$ に含まれる素数をすべて $s \in S$ と $t \in T$ の和で表すことができるような $S, T$ を考ます。これは色々なとり方がありますが、 $d=\lfloor\sqrt{\rm{B2}-\rm{B1}}\rfloor$ を用いて $S=\{\sigma\in [0, d) \ \vert\ \sigma \bot d \}, T=[\lfloor\frac{B1}{d}\rfloor, \lfloor\frac{B2}{d}\rfloor)$ することにします。これはbaby-step giant-step法の考え方と近いです。これを使えばある $\sigma \in S, \tau \in T$ を用いて $pQ=(\sigma+\tau)Q=\sigma Q+\tau Q$ と表せます。なので、 $\sigma Q, \tau Q$ をすべて前計算した後に、 $\gcd((\sigma Q+\tau Q)_Z, N)$ を毎度計算してあげればよいことになります。

ここで、 $\gcd((\sigma Q+\tau Q)_z, N)\neq 0$ は $\sigma Q+\tau Q\equiv 0_{E(\mathbb{Z}/p\mathbb{Z})} \pmod p$ を表していることを思い出します。このとき $\sigma Q\equiv -\tau Q \pmod p$ であるため、Affine座標（射影座標でない、(x, y)の2つ組でとる座標）での楕円曲線の性質を思い出せば $(\sigma Q)_x\equiv (-\tau Q)_x\equiv (\tau Q)_x \pmod p$ （ $P_x$ はアフィン座標でのxをとる記法とする）であることが分かります。まとめると、 $\sigma_x:=(\sigma Q)_x, S_x:=\{\tau_x\ \vert\ \tau \in S \}$ ( $\tau_x, T_x$ についても同様)と表記することとすれば、 $\gcd((\sigma Q+\tau Q)_Z, N)\neq 1$ ならば $\gcd(\sigma_x-\tau_x, N)\neq 1$ であるということです。

また、 $(-\tau Q)_x\equiv (\tau Q)_x$ であるという性質を考えます。すると、 $\sigma Q+\tau Q$ について調べた場合、自動的に $\sigma Q-\tau Q$ についても調べたことになると分かります。これを用いると、アルゴリズムでカバーする集合のサイズはそのままに、片方の集合のサイズを半分に減らすことができます。

この高速化では計算量は定数倍ほどしか変化していませんが、問題を扱いやすい形に整えることができました。ここまでを実装すると、以下のとおりとなります。

def ecm_stage2(n, E, Q, B1, B2):
  d = floor(sqrt(B2 - B1))
  
  S = list([s for s in range((d+1)//2) if gcd(s, d) == 1])
  T = list(range(B1//d*d, B2, d))
  sxs = [(s*Q).to_affine().x for s in S] # x coords of sQ
  txs = [(t*Q).to_affine().x for t in T] # x coords of tQ
 
  for sx in sxs:  
    for tx in txs:    
      if gcd(sx - tx, n) != 1:
        return gcd(sx - tx, n)
  return None

高速化手法2-3: FFT continuation

以上のMontgomery Continuationは計算量が $O(d^2 n\log n)$ でしたが、これを高速化します。 $\sigma_x \in S_x$ を解として持つ多項式 $F$ を考えます。このとき、 $\tau_x \in T_x$ のうち、どれかは $\gcd(F(\tau_x), n)\neq 1$ を満たします。すべての $\tau_x \in T_x$ について $F(\tau_x)$ を求めることは、多項式の多点評価を行うことで $O(d\log^2 d)$ の計算量で行うことができます。 $F$ も分割統治的に掛け合わせること $O(d\log^2 d)$ の計算量で行えるため、この計算量で求めることができると分かりました。

また、 $\tau_x \in T_x$ を解として持つ多項式 $G$ を考えた上で $F, G$ の多項式GCDを取る手法もあります。原理は理解できていないですが、多項式 $F, G$ の多項式GCDを取ると $n$ と互いに素でない値が定数項に出てくるようです。直感的には $\bmod p$ で同じ解を持つことが影響してきそうです。

def polyeval_method(sxs, txs, n):
  F = prod([(sx - X) for sx in sxs])
  for Ftx in multipoint_evaluation(F, txs): # equivalent to [F(tx) for tx in txs]
    p = gcd(prod(Ftxs), n)
    if p != 1: return p
  return 1

def polygcd_method(sxs, txs, n):
  F = prod([(sx - X) for sx in sxs]) # sxs
  G = prod([(tx - X) for tx in txs]) # sxs
  return gcd(polygcd(F, G)[0], n)

高速化手法3: Brent-Suyama’s extension

今までのstage 2では、 $\sigma+\tau$ を掛け合わせていました。Brent-Suyama’s extensionは、掛け合わせる数 $\sigma+\tau$ をある多項式 $f$ を用いて $f(\sigma)+f(\tau)$ とする手法です。 $f$ は低次の $x^e$ やDickson polynomialとされることが多くあります。このように $f$ をとった $f(\sigma)+f(\tau)$ は $\sigma+\tau$ は割り切れる上に他の素因数も含みます。よって、アルゴリズムの成功確率が少し上昇します。

高速化手法4: Edwards Curve

ここからの手法は、20 years of ECMにも記載がなく、gmp-ecmに実装されていない手法です。

今までの手法では、ランダムな楕円曲線を $y^2=x^3+ax+b$ のパラメータ $a, b$ をランダムに決めることでとってきていました。この選択する楕円曲線を、計算するのに都合のよいものに変更することはできないのでしょうか？ECM using Edwards curvesでは、選択する楕円曲線をEdwards Curveと呼ばれる楕円曲線にしました。Edwards Curveでは、射影座標での加算のような制約を受けずにある程度高速な加算が可能な楕円曲線です*8。

また、最新の研究では複数の楕円曲線の表現間を移行しながら演算することで、複数の表現のいいところ取りをするといった手法も提案されています。
eprint.iacr.org

高速化手法5: stage 1 batching

いままでのstage 1では、各素数ごとに最適な加法連鎖を得て、それを用いて乗算を行っていました。ECM using Edwards curvesでは複数の素数をかけ合わせてある程度大きな数にした後に、加法連鎖を計算するといった最適化も加わっています。

以下のでは、掛け合わせる素数の組を工夫し、より演算回数が少ないような組み合わせを探索しています。

eprint.iacr.org

例えば、bound1=256とした際には、以下のように組み合わせることでより効率的に計算することができます(表はFaster cofactorization with ECM using mixed representationsより引用)。

おわりに

時間がなくて最後のほうが雑になってしまいましたが、時間があるときを見て追記していけたらよいかなと思っています。

明日はarkさんによる、about:blankテクニックと"disconnection"の解説です。お楽しみに〜

*1:合成数でも話の大筋は異なりません

*2:楕円曲線は体に対して定義されるものなので、「ようなもの」と表記しています。「ようなもの」なので、加法に対して閉じていないため群をなしません。

*3:循環論法になっていますが、アルゴリズムの正当性を証明することを目指したものではないので目を瞑ってください。

*4:スカラー倍は加算の繰り返しで計算されることから直ちに分かるかと思います。

*5:この分布は $\sin^2 \theta(0\le \theta \le \pi)$ の形をとるようです(Sato-Tate conjecture / 実験コード)。私は一様だと勘違いしており、kobae964さんにご指摘いただきました。ありがとうございます。

*6:これがpower-smoothについての議論でなくて良い理由はあまりよく分かっていません。どこを見てもsmoothnessで議論していて謎です

*7:私の考えでは、Lucas Chainが既知でない場合は $p^k$ に対してPRACアルゴリズムを適用し、それに沿って計算するほうがよいように思います。しかし、gmp-ecmの実装(ecm.c#L1207-1240)ではどうやらそうはなっていません。なぜなっていないのかの理由は謎なので、分かったら教えてほしいです。

*8:Ed25519のEdもこのタイプの楕円曲線を用いていることから名前が取られています