楕円曲線を用いた素因数分解法
はじめに
これは 物工/計数 Advent Calendar 2022 の 6 日目の記事です．
はじめまして．計数 3 年の佐藤 ( @u_zennei ) です．この記事では，RSA 暗号に対する攻撃手法として，楕円曲線素因数分解法 (および，その進化形) を紹介します．
あらかじめことわっておくと，攻撃手法と言っても，* 実際に RSA 暗号を破るような危険性はありません*．実際に使われている秘密鍵がこの手法によって攻撃されうる確率は無視できるほど小さいので，安心して読んでください．
この記事を書いた目的は，攻撃手法を紹介すること自体よりも，代数学が実際にどのように応用されているかを感じてもらうことにあります (そのため，代数数理工学の内容をベースに読めるように書きました) ．
なお，私は未だ学部生として修行中の身であり，記事中にも至らぬ点が多々あるかと思います．ご指摘やご意見をドシドシ頂ければ幸いです．
RSA 暗号と素因数分解アルゴリズム
RSA 暗号は，大きな整数の素因数分解が困難であることを利用した公開鍵暗号方式です．大きな整数が高速に素因数分解できれば，RSA 暗号は破れることになります．
……というだけではあまりにも味気ないので，やや駆け足ですが RSA 暗号と素因数分解アルゴリズムについて紹介したいと思います．
RSA 暗号のしくみ
RSA 暗号は公開鍵暗号の一種で，大きな素数同士をかけあわせることは簡単だが，素因数分解は難しい，という (経験的) 事実を利用しています．以下に，その手順を述べます．
Alice が Bob へ秘密のメッセージを送ろうとしている状況を考えます．Bob は，あらかじめ以下の手順によって自身の公開鍵と秘密鍵を生成しておきます．
- 十分大きな素数 $p, q$ をランダムに選び，$n = pq$ とする．
- 自然数 $e$ を，$1 < e < (p - 1)(q - 1), \; \mathrm{gcd}(e, (p - 1)(q - 1)) = 1$ を満たすように取る．
- 自然数 $d$ を，$1 < d < (p - 1)(q - 1), \; de = 1 \bmod (p - 1)(q - 1)$ を満たすように取る．
- $d$ を自身の秘密鍵として保存する．$(n, e)$ を自身の公開鍵として，Alice に渡す．
Alice は，Bob へ送りたいメッセージを，$n$ 未満の自然数 $m$ に (何らかのきまりで) エンコードします．そして，$c = m^e \bmod n$ によって暗号文 $c$ を計算し，Bob へ送信します．Bob は，$m = c^d \bmod n$ によって平文 $m$ を復元できます．
何気なく書いていますが，この手順が可能であることにも代数の知識を結構使います．
- まず，$(p - 1)(q - 1)$ とは乗法群 $(\mathbb{Z}/n\mathbb{Z})^*$ の群位数です．$\mathrm{gcd}(e, (p - 1)(q - 1)) = 1$ という制約は $e \in (\mathbb{Z} / \varphi(n) \mathbb{Z})^*$ の元であるために必要です ($e < \varphi(n)$ でなくてもよいのですが，計算量が少なくなるように小さく取ります)．
- そして，$d$ は Bezout の公式によって $\varphi(n)$ 未満で一意に定まります．それを求めるには拡張 Euclid 互除法を使って $de + \varphi(n)k = 1$ を解けばよいです．
- 最後に，Bob が $m = c^d \bmod n$ により平文を復元できることは，$c^d = m^{ed} = m^{1 + k\varphi (n)} = m \bmod n$ が中国剰余定理と Fermat の小定理により成り立つことを用いています．
さて，悪意ある第三者が通信路上で暗号文 $c$ を入手したとき，平文 $m$ を復元するにはどうすればよいでしょうか？ パッと思いつくのは，
- $c = m^e \bmod n$ を満たす $m$ を，頑張って見つける ($e$ は公開鍵なので，第三者も知っているとしてよい)．
- 秘密鍵 $d$ を頑張って見つけて，Bob と同じように $m = c^d \bmod n$ として求める．
の 2 つです．特に，$n$ の素因数 $p, q$ を多項式時間で計算できれば，$(p - 1)(q - 1)$ も直ちに求まるので，$de = 1 \bmod (p - 1)(q - 1)$ を拡張 Euclid 互除法により解くことで $d$ が見つかります．すなわち，高速に自然数を素因数分解できるアルゴリズムさえ見つかれば，RSA 暗号は本質的に破られるのです．「RSA 暗号は素因数分解の困難さを利用している」と言われる所以はここにあります．
素因数分解アルゴリズム
現時点では，$c = m^e \bmod n$ から $m$ を求めることも，秘密鍵 $d$ を解読することも，多項式時間で行えるアルゴリズムは見つかっていません．特に，素因数分解の多項式時間アルゴリズムは見つかっていません．
(なお，単純に $2, 3, 5, 7, \ldots $ と割っていけば $O(\sqrt{n})$ で素因数が見つかる (試し割り法) から，多項式時間では？ という気もしてきますが，違います．計算量は入力のサイズに対する関数で考えます．自然数 $n$ を素因数分解するときの入力サイズは $\log n$ なので，試し割り法は指数時間アルゴリズムです．)
素因数分解さえできれば RSA は事実上破れるので，現在に至るまで，高速な素因数分解アルゴリズムが探究されてきました．現状で最速の素因数分解アルゴリズムは一般数体ふるい法で，準指数時間 $\mathrm{O}(\exp( c (\log n)^{1/3} (\log \log n)^{2/3 } ))$ で動きます．その次に，複数多項式二次ふるい法，楕円曲線法と続きます．
楕円曲線法は現時点で 3 番目に速いアルゴリズムで，準指数時間 $\mathrm{O} (\exp (c (\log p)^{1/2} (\log \log p)^{1/2}))$ で動きます．数体ふるい法・複数多項式二次ふるい法と異なり，計算時間が $n$ ではなく $p$ に依存することが特徴です．もっとも，RSA 暗号では素因数 $p, q$ も十分大きいので，この特徴はあまり活かされないのですが．
以下では，楕円曲線法について詳しく述べます．
楕円曲線素因数分解法 (ECM)
楕円曲線素因数分解法 (まどろっこしいので，以下 Elliptic Curve Method の頭字語である ECM と呼びます) は，ごく大雑把にいうと，楕円曲線の Abel 群に対する Lagrange の定理を用いて素因数分解を行なうアルゴリズムです．ECM は H. W. Lenstra によって提案されたもので，古典的なアルゴリズムの部類に入ると思います (完全に余談ですが，H. W. Lenstra には 3 人の兄弟がいて，その全員が数学や数理工学の研究者です)．
ECM について説明する前に，前提となる楕円曲線の概念を導入します．楕円曲線とは大雑把に言えばある種の平面曲線で，曲線上の点同士の「加法」を定めることで Abel 群の構造が入ります．
楕円曲線の Abel 群
楕円曲線は射影平面上で定義されるので，まずその辺りの道具立てを定義します．射影平面とは，直感的には，空間において原点を通る直線上の点を全て同一視して得られる平面のことです．
(このあたりの定義は最終的には使わないので，飛ばして定義 5 以降を読まれても差し支えないと思います)
$K$ を体として，$K^{n + 1}$ なる空間を考える．$P, Q \in K^{n + 1}, \; \lambda (\neq 0) \in K, \; P = \lambda Q$ のとき，$P \sim Q$ として同値関係を定める．$\sim$ による同値類全体がなす集合を $n$ 次元射影空間と呼ぶ．特に $n = 2$ のとき，射影平面と呼ぶ．
射影平面の各点は，アフィン平面 (通常の $K^2$ の点全体) に「無限遠点」を加えたものと自然に同一視できます．
$\alpha, \beta, \gamma \in K, \; \gamma \neq 0$ とする．$(\alpha, \beta, \gamma) \in K^3$ の $\sim$ による同値類を $[\alpha; \beta; \gamma]$ と書き，これを射影座標と呼ぶ．このとき，$[\alpha; \beta; \gamma] = [\alpha / \gamma; \beta / \gamma; 1]$ であり，この点は $K^2$ の点 $(\alpha / \gamma, \beta / \gamma)$ と同一視できる．また，$[\alpha, \beta, 0]$ は無限遠点に対応する．
さて，この射影座標を用いて，楕円曲線を定義します．
$K$ を体，$a_1, a_2, a_3, a_4, a_6 \in K$ とする．$K$ 上の楕円曲線 $E/K$ とは，射影座標 $[X; Y; Z]$ を用いた方程式
$$
X^3 + a_2X^2Z + a_4XZ^2 + a_6Z^3 − Y^2Z − a_1XYZ − a_3YZ^2 = 0
$$
を満たす点の全体である．
この定義は，アフィン座標を用いて書き直すこともできます．そのためには，「分母を払う」処理をしたあとで，無限遠点を付け加えればよいです．
$K$ を体，$a_1, a_2, a_3, a_4, a_6 \in K$ とする．$K$ 上の楕円曲線 $E/K$ とは，アフィン座標 $(x, y)$ を用いた方程式
$$
x^3 + a_2x^2 + a_4x + a_6 − y^2 − a_1xy − a_3y = 0
$$
を満たす点の全体に，** 無限遠点 $\mathcal{O}$ ** を付け加えたものである．
更に，体 $K$ の標数が $2, 3$ ではないとき，更にきれいな形で楕円曲線を表せます．
$K$ を標数が $2$ でも $3$ でもない体，$a, b \in K, \; 4a^3 + 27b^2 \neq 0$ とする．$K$ 上の楕円曲線 $E/K$ とは，アフィン座標 $(x, y)$ を用いた方程式
$$
y^2 = x^3 + ax + b
$$
を満たす点の全体に，** 無限遠点 $\mathcal{O}$ ** を付け加えたものである．
この方程式による定義が，おそらく楕円曲線の定義として最もポピュラーだと思います．せっかく射影平面を定義しておいて申し訳ないですが，以降は楕円曲線を上記のアフィン座標による Weierstrass 標準形で扱うことにします．主に有限体 $\mathbb{F}_p$ (ただし，$p \neq 2, 3$) 上の楕円曲線を扱うので，一つ例を出しておきましょう．
$\mathbb{F}_5$ 上の楕円曲線 $E_1 / \mathbb{F}_5$ を Weierstrass 標準形
$$
y^2 = x^3 + 2x + 1
$$
により定義すると，
$$
E_1 / \mathbb{F}_5 = \{ (0, 1), (0, 4), (1, 2), (1, 3), (3, 2), (3, 3), \mathcal{O} \}.
$$
さて，楕円曲線には，点同士の加法が入ります．加法の定義は多分にテクニカルです．$K = \mathbb{R}$ のとき，$P, Q \in E / K$ の和 $P + Q$ を，「$P$ と $Q$ を通る直線が楕円曲線と交わる点を，$x$ 軸に対して折り返したもの」で定めます．$P = Q$ のときは，$P$ を通る接線を代わりに考えます．この操作を，実数体に限らない任意の体について一般化して加法を定めます．
$K$ を標数が $2$ でも $3$ でもない体とし，$E / K$ を $K$ 上の楕円曲線とする．$P = (x_1, y_1), Q = (x_2, y_2) \in E / K$ の和を以下のように定める．
- $P \neq Q$ のとき．$\lambda = (y_2 - y_1) / (x_2 - x_1), \; \nu = y_1 - \lambda x_1$ として，$P + Q = (\lambda^2 - x_1 - x_2, -\lambda x_3 - \nu)$ とする．ただし，$\lambda$ の分母が $0$ になるなら，$P + Q = \mathcal{O}$ とする．
- $P = Q$ のとき．$\lambda = (3 x_1^2 + a) / 2y_1, \; \nu = y_1 - \lambda x_1$ として，$P + Q = 2P = (\lambda^2 - 2x_1, -\lambda x_3 - \nu)$ とする．ただし，$\lambda$ の分母が $0$ になるなら，$2P = \mathcal{O}$ とする．
この加法により，楕円曲線 $E / K$ に Abel 群 の構造が入ります．単位元は $\mathcal{O}$ です．実は，きちんと Abel 群になることの証明はかなり大変 (例えば，(縫田, 2020) ではこの証明に 13 ページを割いています) なので，ここでは省略します．
こうして楕円曲線の Abel 群が得られたわけですが，群構造はかなり非自明です．例えば，Weierstrass 標準形の方程式 $y^2 = x^3 + ax + b$ が与えられたときに，有理点 (楕円曲線上の点のことをこう呼びます．$\mathbb{Q}$ 上の楕円曲線がよく研究されていたことの名残りだと思います) を全列挙することは相当難しいです．$\mathbb{Q}$ 上では有理点がそもそも存在するのか，存在するとして有限個か無限個か (つまり，群が無限か有限か)，ということを判定するのですら難しいようです．一方で，$\mathbb{F}_p$ 上の楕円曲線については，少なくとも有理点が存在することは保証されます．しかし有理点を具体的に見つけ出すことは難しいです．群位数も explicit には求まらないのですが，幸いにも群位数を多項式時間で計算するアルゴリズムは存在します (Schoof-Elkies-Atkin のアルゴリズム)．ともあれ，こうした構造の複雑さが暗号では盛んに利用されます．
以下で述べる ECM の他にも，楕円曲線暗号や楕円曲線 Diffie-Hellman 鍵交換 (ECDHE) に楕円曲線が利用されています．これらは，楕円曲線上の離散対数問題 (秘密の自然数 $n$ があり，有理点 $P$ と $nP$ だけ与えられたときに，$n$ を復元する問題) が難しいことを利用しています．最近では，耐量子計算機暗号 (e.g. CSIDH) において楕円曲線が使われているようです (耐量子計算機暗号については，(縫田, 2020) が参考になります)．
さて，ECM は，有限体上の楕円曲線の 群位数 がある限られた範囲で (ほとんど) 一様に分布することを用いています．その裏付けとなる定理を (Silverman, 2009) より紹介します．
$p$ を素数，$E / \mathbb{F}_p$ を楕円曲線とする．このとき，
$$
|\# E / \mathbb{F}_p - (p + 1)| \le 2\sqrt{p}.
$$
ただし，$\# E / \mathbb{F}_p$ でその群位数を表す．
つまり，有限体上の楕円曲線の点はだいたい $p + 1$ 個存在し，「誤差」は $2\sqrt{p}$ を超えないということです．更に，位数は $p + 1 \pm \sqrt{p}$ の範囲においてほぼ一様に分布することが知られています．
ECM の流れ
ECM は，上述の位数に関する事実と，Lagrange の定理を組み合わせて利用しています．この節では，アルゴリズムの流れと，それが動作する理由を説明します．
以下，$p, q$ を相異なる素数，$n = pq$ とします．$n$ の非自明な約数 (すなわち $p, q$) さえ見つければ，再帰的に素因数分解が可能です．そこで，素因数分解を直接行なう代わりに，入力された $n$ の非自明な約数を見つけることを考えます．
ECM の処理の流れは以下のようになります (擬似コード環境がなさそうなので箇条書きにしますが，気持ちは擬似コードのつもりです)．入力を $n$ とします．
- $L$ を十分大きな自然数とし，$k = \mathrm{LCM} (1, 2, \ldots , L)$ とする．
- $\mathbb{Z} / n\mathbb{Z}$ 上の楕円曲線 $E / (\mathbb{Z} / n\mathbb{Z})$ と，その上の点 $P$ をランダムに選ぶ ($P$ を初期点と呼ぶ)．
- $kP$ を計算する．その過程で非自明な約数が見つかることがある．見つかればそこで終了し，見つからなければ $E / (\mathbb{Z} / n \mathbb{Z})$ と $P$ を選びなおして繰り返す．
たった 3 項の箇条書きですが，怪しい点がいくつかあります．
まず，$\mathbb{Z} / n \mathbb{Z}$ 上の楕円曲線というのは，これまでの定義から外れています．$n$ は合成数であり，$\mathbb{Z} / n \mathbb{Z}$ は体にならないからです (楕円曲線は体上で定義されていたことを思い出してください)．しかし，ここはグッと我慢して，そのようなセッティングが可能だと認めてください (環上での楕円曲線の定義を厳密に正当化するにはスキーム論を使うらしく，私は不勉強にして詳細を把握していません)．楕円曲線を Weierstrass 標準形で定義すれば，有理点は方程式を$\bmod N$ で満たす座標だと思えばよく，点の加法も $\mathbb{F}_p$ におけるのと同じように$\bmod N$ で考えればよいです．唯一問題が生じるのは，点の加法における $\mathbb{Z} / n\mathbb{Z}$ での除法の存在ですが，これについては後述します．
次に，楕円曲線上の初期点 $P$ をランダムに選ぶという処理はそう簡単ではありません (先述の通り，与えられた楕円曲線の有理点を具体的に見つけることは難しいため)．そこで，実際には，先に点 $P$ の座標 $(x, y) \in (\mathbb{Z} / n \mathbb{Z})^2$ を決めておき，$a \in \mathbb{Z} / n\mathbb{Z}$ をランダムに選んだあとで，$b = y^2 - x^3 - ax \bmod N$ として楕円曲線 $y^2 = x^3 + ax + b$ を取ることにします．すなわち，初期点の座標を決め打ちしておいて，それがきちんと有理点になるような楕円曲線を (後出しで) 構成するわけです．
また，$kP$ の計算ですが，真面目に $P + P$ を $k$ 回繰り返すのは非効率なので，繰り返し二乗法 (楕円曲線については double-and-add method という呼ばれ方が多いです) を利用します．$i = 1, \ldots , \lfloor \lg k \rfloor$ に対して $2^i P$ を計算しておき，$kP$ を復元します．
最後に，「$kP$ を計算する過程で $n$ の非自明な約数が見つかる」という点について説明します (ここが ECM の本質です)．Lagrange の定理により，任意の初期点 $P$ に対して，
$$
\# (E / (\mathbb{Z} / n\mathbb{Z})) P = \mathcal{O}
$$
が成り立ちます (「群位数」倍すると単位元になる)．そのため，$k$ が $\# (E / (\mathbb{Z} / n\mathbb{Z}))$ の倍数であれば，$kP = \mathcal{O}$ となります．
このとき，点の計算において何が起きるでしょうか？ 無限遠点が得られるということは，座標計算のどこかで「ゼロ割り」が生じているはずです．点の和を計算するとき，$\mathbb{Z} / n\mathbb{Z}$ 上の割り算が必要なのは，$\lambda = (y_2 - y_1) / (x_2 - x_1)$ の計算においてのみです．ゼロ割りが生じるのは，$\mathrm{gcd} (x_2 - x_1, n) \neq 1$ のとき，かつそのときに限られます ($\mathbb{Z} / n\mathbb{Z}$ における可逆元は $n$ と互いに素な数であることを思い出してください)．
さて，こうなるともはや $kP$ の計算は進められない代わりに，$\mathrm{gcd} (x_2 - x_1, n) \neq 1$ を満たす整数 $x_2 - x_1$ が得られています．このとき，$\mathrm{gcd} (x_2 - x_1, n) = p, q, n$ のいずれかですが，$n$ となるのはよほど運の悪い場合に限られます．したがって，$\mathrm{gcd}$ を計算することで $p$ か $q$ が得られます．これで ECM が成功します．
以上の説明に基づくと，ランダムに選んだ楕円曲線の位数が $k$ を割り切る (Lagrange の定理が使える) ことこそ，ECM が成功する条件であることが分かります．$k$ は $1, \ldots , L$ を約数として持つので，楕円曲線の位数が十分小さい素数同士の積で書けるならば，その条件が満たされます．ここで，楕円曲線の位数がある範囲で一様に分布することが効いてきます．もし群位数が特定の値，例えば $p + 1$ しか取らなかったならば，運悪く $p + 1$ が非常に大きな素因数を含んでいたとき，$k$ をものすごく大きく取らなければなりません．しかし，実際には群位数は十分ランダムなので，何度か楕円曲線を取り替えていれば，いつかは群位数が $k$ を割り切るような曲線を見つけられます．これが，ECM が素因数を見つけられる理屈です．
実装と例
拙いですが，ECM を C++ により実装したものを GitHub に置いています．これを使って数値例を示します．
ECM により，$n = 2183 = 59 \times 37$ の素因数分解を試してみます．$L = 5, k = 60$ とします．
まず，初期点 $P = (1982, 1507)$ を取り，楕円曲線 $E/(\mathbb{Z} / 2183 \mathbb{Z})$ を $y^2 = x^3 + 1461 x - 474$ により定めます．$kP$ を計算するために，$2P, 2^2 P, \ldots , 2^5 P$ までを計算すると，それぞれの座標は以下のようになります．
| $2P$ | $(105, 1106)$ |
| $2^2 P$ | $(1377, 280)$ |
| $2^3 P$ | $(1933, 492)$ |
| $2^4 P$ | $(1784, 534)$ |
| $2^5 P$ | $(379, 1580)$ |
これらから，$60P = 2^2 P + 2^3 P + 2^4 P + 2^5 P$ を計算すると，
| $2^2P + 2^3P$ | $(1999, 1490)$ |
| 上$ + 2^4 P$ | $(815, 471)$ |
| 上$ + 2^5 P$ | $(960, 851)$ |
となって，$kP$ が求まってしまいました．曲線を取り直してやりなおします．初期点 $P = (1996, 2149)$ を取り，楕円曲線 $E/(\mathbb{Z} / 2183 \mathbb{Z})$ を $y^2 = x^3 + 343x - 1258$ により定めます．このとき，
| $2P$ | $(1148, 1799)$ |
| $2^2 P$ | $(408, 1198)$ |
| $2^3 P$ | $(1148, 2095)$ |
| $2^4 P$ | $(408, 1457)$ |
| $2^5 P$ | $(1148, 1799)$ |
と求まりますが，$60P$ を計算するために $2^2 P + 2^3 P$ の座標を考えると，$\lambda$ の分母 $x_2 - x_1 = 740$ であり，$\mathrm{gcd}(740, 2183) = 37$ となって，$n$ の素因数が見つかりました．この曲線では $kP = \mathcal{O}$ が成立したことになります．
ECM の計算時間について，以前手元で簡単に試した結果を述べておきます．AMD Ryzen 7 1700 $8$ コア $16$ スレッドの PC を $7$ つ束ね，$56$ ノードの並列計算環境で ECM を走らせたところ，$40$ 桁 ($133$ bit) の合成数から素因数を見つけるのに $2.1$ 秒，$60$ 桁 ($200$ bit) の合成数から素因数を見つけるのに約 $165000$ 秒 (約 46 時間) かかりました．
これは素朴なアルゴリズムと比べると桁違いに速い (例えば，試し割り法で $60$ 桁の合成数を素因数分解するのは，一生かかっても無理でしょう) のですが，それでも RSA 暗号の解読は絶望的です．現在は RSA 暗号で使う合成数のサイズとして $2048$ bit 以上が推奨されています．よほど脆弱な鍵でも $1024$ bit 以上が使われています．$200$ bit の素因数分解に $2$ 日かかるようでは，$1024$ bit の素因数分解は無理でしょう．
ECM より速い 2 つのアルゴリズム (一般数体ふるい法，複数多項式二次ふるい法) では，これよりも計算時間が改善しますが，それでも本質的には解読は難しいと考えられています．現在までに素因数分解された最大の合成数は $829$ bitで，要した計算時間は $2,700$ core-years だそうです (Boudot et al., 2020)．今のところ我々は，鍵に使う合成数を十分大きく取れば，RSA 暗号を安全に使うことができます．(大昔に生成した鍵長 $800$ bit 以下の鍵を使いつづけている人は今すぐ更新するべきだ，ということでもあります；ssh-keygen は $1024$ bit 未満の鍵を生成できないので，よほどレアなケースにはなりますが……)
一方で，極めて特殊な条件下では，今まで述べてきた一般的なアルゴリズムとは別の手法により，現実的な時間での素因数分解が可能です．例えば，素因数 $p$ と $q$ が非常に近い場合などは，「特殊な条件」の一つです．このとき，$\sqrt{n}$ の十分近くに限って試し割りを繰り返せば素因数が見つかってしまいます．以下では，「特殊な条件」に限った攻撃手法の一つとして，素因数 $p$ が $4p = 1 + Dv^2$ と書ける場合に対する素因数分解法を述べます．
虚数乗法を用いた ECM
理論
ECM のポイントは，「位数が小さな素数の小さな冪の積となっている楕円曲線があれば，その上で点を $k$ 倍することで，Lagrange の定理によりゼロ割りが発生し，素因数分解ができる」ことでした．
さて，ここでもし群位数が既知の楕円曲線 $E / (\mathbb{Z} / n\mathbb{Z})$ があったらどうでしょうか？ もしそんな楕円曲線があれば，適当な点を取って，ちょうど群位数の回数だけ足し上げることで，ただちに $n$ の素因数が求まるでしょう．楕円曲線を何度も取りなおして，運良く群位数が $k$ を割り切るようになるまで試す必要がなくなります．そして，点の加法は double-and-add method により群位数の対数オーダーで実行できますし，群位数は Hasse の定理によりたかだか $p + 1 + 2\sqrt{p}$ なので，全体として $p$ の多項式時間で素因数分解が走ります．
そして，$p$ が特定の条件を満たすなら，群位数が既知の楕円曲線を実際に構成できます．その条件とは，$D \in \{3, 4, 7, 8, 11, 19, 43, 67, 163\}, \; v \in \mathbb{N}$ により，
$$
4p = 1 + Dv^2
$$
と書けることです．このとき，位数がちょうど $p$ の楕円曲線 $E / (\mathbb{Z} / n \mathbb{Z})$ を $1/2$ の確率で構成することができます (具体的には，Hilbert 類多項式 $H_D (x)$ の根を $j$-invariant に持つ楕円曲線を取ればよい)．
この記事の範囲で理屈を説明するのはとても無理なので，気持ちだけサラッと触れておきます (ここは読み飛ばしてもらっても大丈夫です)．上述の条件が満たされるとき，
$$
p = \frac{1 + \sqrt{-D} v}{2} \frac{1 - \sqrt{-D} v}{2}
$$
と書けます．このとき，虚二次体 $\mathbb{Q} (\sqrt{-D})$ のある部分環による虚数乗法を持つ (自己準同型環が $\mathbb{Q} (\sqrt{-D})$ のある部分環と同型になる) 楕円曲線が存在します．さて，$t = \# E - (p + 1)$ は Frobenius 自己準同型と呼ばれる写像のトレースですが，$\mathbb{Q} (\sqrt{-D})$ において $(1 \pm \sqrt{-D} v) / 2$ はノルムが $p$ でトレースが $\pm 1$ の元となっています．このとき，$t = \pm 1$ が $1/2$ ずつの確率で成り立ちます．$t = 1$ であれば，$\# E = p + 1 - t = p$ が従います．
(なお，$D$ のリストがマジックナンバーすぎると思った人もいるかもしれませんが，これらの数 ($3, 4, 7, \ldots , 163$) には Heegner 数という名前がついています．詳しくは「虚二次体の類数」といったワードで調べてみてください．)
さて，位数 $p$ の楕円曲線が得られたので，あとは適当な初期点 $P$ を取って $n$ 倍すれば，ECM の原理によって素因数分解に成功します．問題は有理点 $P$ を取るのが難しいことです．従来の ECM では，初期点を先に決めて，楕円曲線を「後出し」することで問題を回避していました．しかし，使う楕円曲線が先に決まっている以上，この手は使えません．
そこで，以下のような手法を取ります ( (白勢, 2016) による )．得られた位数 $p$ の楕円曲線の方程式を $y^2 = x^3 + ax + b$ とします．初期点 $P$ の $x$ 座標を適当に決め打ちし，$x_0$ とします．$\tau = x_0^3 + ax_0 + b$ が $\mathbb{F}_p$ において平方元である確率は $1/2$ です ($\mathbb{F}_p$ のちょうど半分は平方元)．$p$ が既知ならば $\mathbb{F}_p$ 上で平方根を求めるアルゴリズムによって $y = \sqrt{\tau}$ が見つかるのですが，そうではありません．ここで，多項式環の剰余環
$$
R = (\mathbb{Z} / n\mathbb{Z} [X]) / (X^2 - \tau)
$$
を考えて，その上の楕円曲線 $E/R$ に話を移します．$P = (x_0 + 0X, 0 + X) = (x_0, X)$ は $E/R$ の有理点となっています．そこで $nP$ を計算すると，分母に $p$ の倍数に対応する $R$ の元が生じる (ゼロ除算) ことで素因数が見つかります．「$p$ の倍数に対応する $R$ の元」とはなんでしょうか？ そのような元の特徴付けを試みます．
写像 $\varphi : \mathbb{F}_p [X] / (X^2 - \tau) \to \mathbb{F}_p$ を，$\varphi (f(X)) = f(\sqrt{\tau})$ として定めます (先述の通り $\sqrt{\tau}$ は計算できないのですが，$\varphi$ は議論のために導入した写像で，実際に計算する必要はありません)．これは代入写像であって，well-defined な全射準同型になっています．よって，準同型定理から，
$$
(\mathbb{F}_p [X] / (X^2 - \tau)) / \ker \varphi \simeq \mathbb{F}_p
$$
が成り立ちます．さて，$a_0 + a_1 X$ が $p$ の倍数 (すなわち $0 \in \mathbb{F}_p$) に対応しているとします．このとき，$\varphi (a_0 + a_1 X) = 0$ です (同型は単位元を単位元に移す)．ゆえに，
$$
a_0^2 - a_1^2 \tau = (a_0 + a_1 \sqrt{\tau})(a_0 - a_1 \sqrt{\tau}) = \varphi (a_0 + a_1 X)(a_0 - a_1 \sqrt{\tau}) = 0.
$$
よって，$\sqrt{\tau}$ を計算することなく，$\mathbb{F}_p$ における単位元 ($p$ の倍数) として $a_0^2 - a_1^2 \tau$ を得ることができます．
(補足 : この計算は本質的に終結式と結びついており，終結式の計算を考えることで $D$ の制約を外せる，という研究があります (相川, 2018) ．)
以上の議論をまとめると，$E / \mathbb{Z}/n\mathbb{Z}$ の代わりに，$E / R$ 上で ECM を走らせ，座標の分母に可逆でない元 $a_0 + a_1 X$ が現れたならば，$\mathrm{gcd}(n, a_0^2 - a_1^2 \tau)$ が $n$ の素因数を与えることになります．このアルゴリズムが成功する確率は，選んだ楕円曲線で $t = 1$ となっている確率 ($1/2$) と，$\tau$ が $\mathbb{F}_p$ の平方元となっている確率 ($1/2$) の積なので，$1/4$ となります．
実装
このアルゴリズムの実装は，先述した リポジトリ の CM ブランチに置いています．
この手法では，一回あたりの試行では $\mathrm{O}(\log n)$ 反復しか回らず，各試行は $1/2$ の確率で成功するので，全体の処理は $n$ が数千 bit であっても数秒で終了します．手元の環境では，$2048$ bit の合成数を $3$ 秒弱で素因数分解できました．
しかしながら，最初に述べた通り，これは RSA 暗号への本質的な脅威となりません．その理由は，$p$ が満たすべき条件が厳しすぎるからです．条件を満たす $p$ がどのくらいあるか，ラフに見積もってみましょう．そのためには，$1 + Dv^2$ がどのくらいの頻度で数直線上に並んでいるか見れば十分です．$D$ の取りうる値は有限なので，動かせるパラメータは $v$ のみです．すると，$1 + Dv^2$ は $2$ 乗のオーダーで値が増えていきます．一方で，素数はもっと密に存在します (素数定理より)．実のところ，ランダムに取った素数 $p$ が $4p = 1 + Dv^2$ と書ける確率は，$\sqrt{D / p}$ 程度であるとされています．$p$ は通常 $2^{1024}$ 以上で取るので，その確率は極めて小さいことが分かります (観測可能な宇宙の粒子数は $2^{265}$ 個くらいらしいです．宇宙から粒子をランダムにピックして特定の一つをたまたまツモる確率よりもなお圧倒的に小さいのです)．
おわりに
ここまで読んでいただき，ありがとうございました．
代数学がどのように使われているか，自分の語れる範囲で少し書いてみました (私の不勉強により釈迦に説法をする羽目になっていないか少し不安ですが……)．なんとなくの雰囲気でも伝わったなら嬉しいです．
最後になりましたが，明日以降の 物工/計数 Advent Calendar 2022 もお楽しみに！
