前 6‒06 p ｜ 
もくじ i ｜ 次 6‒08 n 
楕円曲線を使って、巨大整数に含まれる数十桁の因数を検出できる。計算は、曲線上の勝手な点を選んで整数倍するだけ。
楕円曲線法ステージ2は、モンゴメリー、ブレント、日本のスヤマが開拓した面白い領域だが、一般向けの解説（教科書やウェブの記事）があまりない。出発点となるレンストラーのステージ1と、その先にあるモンゴメリー形式、標準版ステージ2、素数ペアリングについて整理してみた。前提として、中学数学の知識を仮定する。
20世紀後半まで「整数を因数分解するには、試しに割ってみるしかない」と考えられていた。1970年代に英国のジョン・ポラード（
John. M. Pollard）が「p − 1」法を発見して、「他にも方法がある」ことを示した。
合成数 n が与えられ、それを因数分解したいとする。つまり n の約数 g を知りたい。例えば n = 1279037 なら 631 × 2027 と分解され、g = 631。小さい素数で片っ端から割る方法以外で、この値を得たい。ポラードの方法は、次の3ステップから成る。
小さな因数をたくさん含む整数 m を考える。典型的には、適当な自然数 B1 を選んで「1, 2, 3, …, B1 の最小公倍数」を m とする。例えば B1 = 10 とすると:
小さな自然数 J を選んで、K ≡ Jm (mod n) を計算する。例えば J = 3 とすると、32520 を n = 1279037 で割った余りが Kとなる。コンピューターを使えば、この計算は簡単にできる:
K − 1 と n の最大公約数を g とする。ユークリッドの互除法（
）を使えば、n の因数分解を使わずに最大公約数を計算できる。
…魔法のように g が求まった。
今考えれば、上記は「公開鍵暗号などで、よくあるタイプの計算」だが、1970年代まで誰もこれを考えなかった。フェルマーやガウスですら、このシンプルな計算を思い付かなかった。
「人類は1970年代に突然頭が良くなった」というわけでは ないだろう。コンピューターがあれば Jm (mod n) は何でもない計算だが、コンピューターがない時代、この種の計算を湯水のように使うという発想は、非現実だった。計算機の発達は人類の意識・世界観を変化させ、「何が現実で何が現実でないか」が変わった。
Jm (mod n) の値を求めるとき、実際に Jm を計算する必要はなく、通常 m を計算する必要すらない。m は小さな素数（または素数べき）の積なので、まとめて m 乗しなくても、少しずつ「小さな素数」乗してその都度 n で割った余りを考えれば同じことになる。一つ一つの「小さな素数」乗の部分は、「2乗して n で割る」という単純な計算を繰り返すことに帰着される（繰り返し2乗法）。「p − 1」法の実装については別記事「[JS] メルセンヌ数の分類と分解」を参照。
因数分解が成り立つ理由は次の通り。
n が、異なる2素数 p, q の積に分解されるとする（n = pq）。実際に因数分解に成功するまで p と q の値は分からない。それでも未知の p について、フェルマーの小定理（）により Jp − 1 ≡ 1 (mod p) が成り立つ（q についても同様）。
従って、もし m が p − 1 の倍数なら、m = a(p − 1) として:
つまり Jm − 1 は p の倍数。
実際には p の値は未知なので Jm (mod p) を計算できないが、K ≡ Jm (mod n) なら計算できる。mod p で ≡ 1 になる Jm は、mod n では何になるか？
これは中国剰余定理（）のシチュエーションであり、答えは mod q で何になるかに依存するが、一般的に言えば、何になるか分からない。特別な場合を除けば、Jm ≡ 1 (mod n) ではなく、従って Jm − 1 は n の倍数ではない。
要するに Jm − 1 は p の倍数だが、n の倍数ではない。言い換えれば、p は Jm − 1 の約数だが、n は そうではない。一方、仮定により p は n の約数。そうすると、Jm − 1 と n の最大公約数を求めれば、それが p になるはずだ。この場合、Jm − 1 の計算は mod n で行って差し支えない。
…言われてみれば単純な話で、フェルマーがこれに気付いたとしても、原理的には少しもおかしくない。しかし前述のように、当時としては これは人類の「現実」の範囲外だった。
このアルゴリズムは、J の選び方が悪いせいで失敗することもあるが、その場合は単に J の値を変えればいい。失敗の原因は、ほとんどの場合、J ではなく B1 の選択に関係している。
B1 は勝手に選んだ数なので、「もし m が p − 1 の倍数なら」という部分が成り立つ保証はない。最初の例では、因数分解
を実行した。この場合 p = 631 と考えることができるが、そうすると、上記の「もし」は「もし m が 631 − 1 の倍数なら」という意味になる。勝手に選んだ数に基づく数値 m が、たまたま 630 の倍数になっている…というのは、ずいぶん虫のいい仮定で、一般論としては、そんなにうまく いくわけない。
ただし p − 1 が比較的小さい場合、この問題は簡単に克服される。考えている m は「B1 以下の自然数の公倍数」なので、B1 以下の自然数の積として表される数は自動的に m の倍数になるからだ。具体的に:
この数は、10 未満の素数（または素数べき）しか含んでいない。
一方、われわれの m は:
この数は、10以下の全ての自然数の倍数なので、当然 630 = 2 × 9 × 5 × 7 の倍数にもなっていて、「もし」が満たされる。
もちろん p − 1 の分解がもっと大きな因数を含む可能性もあり、そうすると B1 が小さ過ぎて、アルゴリズムは失敗する。例えば、下記の21桁の整数の分解を考える:
このとき、p − 1 の分解は次の通り:
上記は後から分かることで、最初は p が不明なので、都合よく B1 を選択することができない。とりあえず B1 = 500 とすると、それに対応する m は 8 の倍数であり、27 の倍数でもあり、5 の倍数でもあり、67 の倍数でもあるが、2677 の倍数には なっていないので、「もし」は成り立たない。では思い切って B1 = 5000 などと、大きくしたらどうか。その場合、m は 2677 の倍数でもあるので（それどころか 4998 の倍数でも 4999 の倍数でもあるので）、自動的に上記の数の倍数にもなり、「もし」が満たされる。
B1 は大きければいい というものではなく、大き過ぎると、p − 1 ばかりか q − 1 も同様の関係を満たしてしまう。その場合 K − 1 は q の倍数にもなってしまい、K − 1 と n の最大公約数は n 自身になってしまう。これでは目的を達成できない。しかし大き過ぎる場合、通常、単に B1 を小さくすればいい。問題は小さ過ぎるケースだ。
大きな素数 p を考えた場合、それがどんな数であろうと、p − 1 を素因数分解すれば比較的小さな整数の積となる（最悪のケースでも 2 × (p − 1)/2）。B1 を「p − 1 の最大素因数」以上にできれば、「もし」は満たされる。例えば p が18桁の素数だとすると、かなりの確率で、p − 1 は「9桁の素数2個の積」か「それよりもっと小さい素因数の積」に分解されるだろう。その場合、B1 を10億（10桁）まで上げられるなら、因数分解したい整数に含まれる因数 p が検出される。
検出したい素数 p がもっと大きい場合（例えば20桁や30桁の場合）、一般には p − 1 の最大素因数は9桁より大きくなる。B1 が10億のとき、20桁の素因数はある程度検出可能かもしれないが、30桁の素因数はほぼ検出できない。理論的には B1 をもっと大きくすれば対応できる範囲は広がるけれど、B1 を1桁大きくするごとに計算量はどんどん増大するので、現実的には限度がある。
p が20桁だが、うまくいく例:
p − 1 の最大素因数が2797万台なので、B1 を3000万程度にできれば、うまくいく。p は「サニレ・サニニ・イココクロシクロシ・サココロク」という面白い数字の並びを持っている。
一方、同じ20桁の因数でも、次の場合はどうか。
p − 1 の最大素因数が37兆台になってしまった。つまり、20桁の因数「ハロロ・コロニ・ロハコロロニハニ・イハサイコイ」（p = 86656268566282183151）をポラードの方法で検出したい場合、B1 をざっと38兆まで上げない限り、うまくいかない。「10億」くらいならともかく「38兆」は、きつい要求だ。20桁の因数でこれでは、因数が25桁・30桁…となったとき困ることは明らか。何とか B1 を値切れないだろうか？
レンストラー（Lenstra）はオランダの数学者兄弟。楕円曲線法を発見したのは、そのうちヘンドリク・レンストラー（Hendrik Willem Lenstra, Jr.（））。1985年2月に成果が発表され、1987年に論文 [Len87] が出版された。
抽象代数的な議論になって しまうけれど、この節では動作原理を簡単にまとめておく。§3 ではもっと具体的に説明する。
1970年代に発見されたポラードの「p − 1」法は画期的だったが、限界もあった。
レンストラーが目を付けたのは、「ポラードの方法がうまくいくのは、乗法群の位数が p − 1 だからだ。この部分は一般化できる」という点だった。どういうことか というと…
この観点からすると、フェルマーの小定理は、ラグランジュの定理の一例となる。つまり「ある整数を p − 1 乗（または p − 1 の倍数乗）すると ≡ 1 (mod p) になる」という数論的性質は、「ある元 a について、 a * a * a * … のように、群の位数回（あるいは位数の倍数回）群の演算を繰り返すと、結果は単位元になる」という群論的性質の現れと考えることができる。そうすると、
という「p − 1」法は、より一般的に、
と言い換えることができる。
これだけだと「p − 1」を「群の位数」と小難しく言い換えただけで、意味がないようだが…
群の位数は p − 1 とは限らない！
当たり前の事実ではあるが、そこに着目することが決定的な一歩となる。
例えば、位数が p − 2 の別の群がどこかに あったとして、その群において同じことを考えた場合、p − 2 が分解しやすければ因数分解が成り立つことになる。もしそれも駄目でも、位数が違う群をいろいろ試せば、それらの位数のどれか一つくらいは、比較的小さい因数しか持たないはずだ。
このような多様な群の供給源として、楕円曲線の点の群（後述）を使えばいい。
これがレンストラーのアイデアだった。
具体例として、20桁の「ハロロ・コロニ」 p = 86656268566282183151 について、p − 1 の分解から攻めようとすると「38兆」という きつい要求が出てくるのだった。ところが、楕円曲線
を考えたとすると、実はこの曲線上の点の数（群の位数）は p − 1 ではなく p − 453604093 になっている。この位数は
と分解される。従って、ポラードの方法と同様のアルゴリズムによれば、B1 = 11000 程度でうまくいく。「38兆」に比べると、ずいぶんリーズナブルだ。
曲線 y2 ≡ x3 + 11x + 1 (mod p) は、どこから出てきたのか？ これは定数項を 1 に固定して、1次の項の係数を 1, 2, 3, 4, … と変えて順に試したら、たまたま 11 で成功しただけ。この曲線でなければいけない わけではなく、この曲線なら常にうまくいく わけでもない。重要なのは、ある曲線が駄目でも別の曲線があるということ。乗法群は位数 p − 1 に固定されていて交渉の余地がないが、楕円曲線は何種類もあるので、位数の相性が悪ければ曲線そのものを乗り換えることができる。
「p − 1」法の核心は、勝手に選んだ数 J の「べき乗」だった:
つまり J を m 回、掛け合わせた。それに対応する楕円曲線バージョンの計算は「勝手に選んだ曲線上の点 J を m 回、足し合わせること」。「掛け合わせる」と「足し合わせる」で用語は違うが、本質的に同じパターンの計算だ。
以下では「点を足し合わせる」という処理を定義して、上記のアイデアを実際に試してみる。
楕円曲線の「点の和」というのは、内容的には単なる四則演算。実際の計算はコンピューターがやってくれるので、人間にとっての計算難易度は問題ではないけれど、一つ一つの足し算は比較的シンプルな処理だ。
楕円曲線は広い意味での3次曲線で、その標準形は:
典型的には下図のような曲線となる。
図の P や Q は楕円曲線上の点。「点と点を足す」ということは、P + Q を例に取ると、図の R を求めることに当たる。足し合わせたい2個の点を通る直線を考えると、この直線は曲線と別の場所でもう一度交わる。その第3の点を用いて、点の和が定義される。「この第3の点を P + Q と定義する」なら単純だが、ひとひねりあって、実際には「その点を −R として、y-座標の符号を逆にした点を R とする」と約束する。
資料によっては、図の R と −R の名前を逆にしている場合もあるが、名前が逆なだけで、定義の内容に違いはない。
曲線上の任意の点 P = (x1, y1) に対して、(x1, −y1) も曲線上の点であり、それらは x-座標が同じで y-座標の符号だけが異なる。後者は、前者から見た「マイナスの点」 −P となる。
P の座標を (x1, y1)、Q の座標を (x2, y2) とした場合、この2点を通る直線の傾きは:
P = Q の場合、上記のままでは λ（ラムダ）の分母が 0 になってしまい計算できないが、「P が Q に少しずつ近付いて極限で一致した」と考えれば、このときの λ は P での接線の傾きと考えるのが妥当だろう。接線の傾きは:
いずれの場合も、R の座標 (x3, y3) は、λ を使って次のように表現される。
これらの式は、直線と曲線の交点を代数的に求めることで得られる。
「接線の傾き」だけは中学数学の範囲（記事の前提）では説明できず、微分が必要になる。ここでは理屈抜きに、上記の公式を「和の定義」としても構わない。
ただし、ここまでは「ゼロ」が絡まない場合の話だ。
ある点 P について −P は x-軸を挟んで反対側の点だが、
は、どこにあるのだろう？ 当然 P + (−P) は「ゼロ」になってほしいが、「ゼロ」はどの点だろうか？
「ゼロは原点」と言いたいところだが、原点は（一般には）曲線上の点ではないし、和についての幾何学的説明にも反している。仕方がないので、強引に
」から成る点の集合を考え、「 も楕円曲線上の点だが、 には座標がない」とする（ が曲線上のどこにあるのか ということは不問に付す）。そして「この は普通のゼロのように振る舞う」と約束する。すなわち、曲線上の任意の点 P について:
= + P = Pこれで「ある一つの楕円曲線の上の普通の点たち＆ゼロと呼ばれる点」という集合の要素間に「点の和」が定義され、逆元（マイナスの点）も定義される。この演算は、結合法則と交換法則を満たす。 は加法の単位元（中立元）となる。
幾何学的には、 は曲線上において無限の遠方にある点（無限遠点）と解釈され、射影平面上では実体のある「まともな点」になる。
ここまでは、点の座標値が実数であるかのように記述してきた。
実際の楕円曲線法では、点の座標値は「整数のようなもの」となる。具体的には a, b を 0 でない整数、n を 3 より大きい奇数として、
を満たす整数 (x, y) をこの曲線上の点と認定する。
ただし座標値について、n で割った余りが同じなら「同じ」と見なす。例えば n = 7 とすると、x = 1 と x = 8 は「同じ」と見なされる。「ある月の1日と8日は、曜日が同じ」というような発想だ。
という式は、このような意味を持つ（ 整数の合同）。
例えば曲線 y2 ≡ x3 + 5x + 1 (mod 7) 上において、点 P = (1, 0) と点 P′ = (8, 0) は、座標値が「同じ」なので、点として等しい。このことを
と書き表すことにする。
この場合にも、 という特別な点が曲線上にあると考える（ の具体的な座標値は不問に付す）。そうすると、前述の公式をそのまま使って「点の和」を考えることができ、n が素数なら、点と点の和が矛盾なく定義される（n が合成数なら、点と点の和が定義されない こともある）。
これは普通の意味での「曲線」ではなく、幾何学的には「碁盤の上に散らばった碁石」のような姿を持つ（→ 補足1）。
次の楕円曲線 E を考える。
P = (0, 1) は、曲線 E 上の点だ（曲線の定義式に x = 0, y = 1 を代入することによって確認できる）。−P = (0, −1) もこの曲線上にあり、座標値についての規約から −P = (0, 6) とも書ける。
試しに、和の定義に従って Q = P + P を計算してみよう:
よって Q = (4, 5) もまた、この曲線上の点であることが分かる。このことは、曲線の定義式に x = 4, y = 5 を代入することにより直接確認できる（両辺とも ≡ 4 になる）。
λ の計算には割り算が絡んでいるが、今は「整数のようなもの」を考えているので、常に割り算ができる保証がない。上記の例では、たまたま 4/2 となって、普通に計算しても結果は合っているが、分子が分母の倍数でないときは困ってしまう。
実は、この座標計算においては「割り算は、分母の逆数（逆元）を掛けること」と解釈される。ただし「逆数とは、もとの数と掛け合わせて結果が ≡ 1 (mod n) になるような数」とする。この例では、2 × 4 = 8 ≡ 1 (mod 7) なので、2 の逆数は 4 であり、これを掛ければいい。
逆数が存在する場合、拡張互除法（ユークリッドの互除法の拡張版）によってそれを求めることができる（§4）。mod 7 での 2 の逆数を手っ取り早く求めるには、端から順に暗算して「ニニンがシ＝違う。ニサンがロク＝違う。ニシがハチ＝7で割ると余り1＝これだ！」とやればいい。
法が素数なら、0 以外の数は必ず逆数を持ち、座標値の四則演算が問題なく実行される。
同じ曲線 E について、異なる点を足す場合の例として、Q + P = (4, 5) + (0, 1) が R = (4, 2) に等しいことを示す。
Q = P + P なのだから、上記の計算は 3P = P + P + P = (4, 2) と書くこともできる。R = (4, 2) = (4, −5) なので、R = −Q でもある。
さらに、R + P = (4, 2) + (0, 1) が S = (0, 6) に等しいことを示す。
上記の計算は 4P = P + P + P + P = (0, 6) に当たる。(0, 6) は −P に等しいので、もう1回 P を足すと（ についての定義から）結果は になる。
結局、曲線 E は、5個の点 P, Q, R, S, を持つことが分かった。±P, ±Q, と書くこともできる。曲線の定義の式の x と y に 0, 1, 2, …, 6 を入れて、全ての組み合わせの可能性を試すと、 E にこれ以外の点はないことが分かる。つまり、この楕円曲線の点の群の位数（要するに点の数）は 5 であり、点 P を 5 回足し合わせると（加法の）単位元 になる。
これはフェルマーの小定理の楕円曲線版のような現象だが、フェルマーの小定理の場合、背景にある乗法群の位数は p − 1 であり、この回数、演算を繰り返すと（乗法の）単位元 1 が得られるのだった。ところが、今考えた楕円曲線の点の群の場合、p = 7 という素数を法としているが、位数は p − 2 になっている。この群を使って「p − 1」法と同様のことを行うなら、「p − 2」法が成立する！
楕円曲線の点の群の位数は p − 2 に固定されている わけではない。曲線を定義する式を
として、パラメーター a を変えた場合、点の数（位数）はかなりランダムに変動する。「何種類も曲線を試せば、そのうち一つくらいは因数分解に都合がいい位数になるだろう」という期待が生じる。
n = pq という素因数分解を実行したいが、p と q が分からないとする。
「p − 1」法の基本戦略は、p − 1 の倍数になっていると期待される m を構成することだった。p − 1 というのは mod p の乗法群の位数に他ならない。期待通りになった場合:
楕円曲線法の基本戦略は、同様に、「楕円曲線の点の群の位数」の倍数になっていると期待される m を構成すること。期待通りになった場合:
(mod p) となる。 (mod q)。ここで J は勝手に選んだ楕円曲線上の点であり、[m]J は、その点を m 回足し合わせたものとする。（角かっこに深い意味はなく、単に mJ のように書かれることも多い。）
「mod p では結果は だが、mod q では結果は ではない」という上記の状況は、mod pq つまり mod n において、どういう状態なのか？
直観的に、mod p で測って無限の遠方にある点 は、上位互換の mod n でも無限になるはずだ。同様に、mod q で測って有限の場所にある点は、mod n でも有限のはずだ。
この二つのことは矛盾している。mod n の視点では、心の半分では「この和は無限」と考え、心の半分では「この和は有限」と考えるので、解決不可能なジレンマが発生する。すなわち、この和は定義不能・計算不能となるはずだ（→ 補足2）。
「計算不能」とは、具体的にどういうことか？
点の和の計算は座標値の四則演算であり、計算できなくなる ような場所はあまりない。しかし一カ所、デリケートな場所があった。λ の計算で、割り算（逆数の計算）が必要になるのだった。法が合成数の場合、必ずしも逆数は定義されないので、そこで計算できなくなる ことがある。逆に、計算できなくなる としたら、そこしかない。つまり「計算不能」の正体は「逆数が定義されない」だ。
k の逆数 k−1 (mod n) を計算しようとして失敗するとき というのは、k と n が互いに素ではないとき、つまり k と n の最大公約数が 1 ではない ときである。その場合、確かに計算できなくなるが、それは「失敗」だろうか？
そもそもの目的は n の因数（約数）を知ることなのだから、n の約数が出てきたら目的は達成される。エラーの原因のその最大公約数が、探していた因数だ！
この「エラー」はコントロールされたものであり、「m が mod p での曲線の位数の倍数だが、mod q においては そうではない」という状況を作ることにより、人為的に引き起こすことが できる。つまり、曲線の位数が比較的小さな素因数しか持たず、それらと比べて B1 が適度に大きいときに発生する。
楕円曲線の点の和を計算する関数を作り、曲線上の勝手な点を m 回足し合わせて、何が起きるか実際に試してみよう。別記事で作ってあるライブラリを流用して JavaScript で記述する。
点の和の計算で逆数 k−1 (mod n) を使うので、まずは拡張互除法を用意する。
以下の実装は [Ste07] の35–37ページに基づく。while ループを抜けたとき、ローカル変数 a は k と n の最大公約数になっていて、この値が 1 なら r = k−1 (mod n) となる。
function Inv_mod( k, n ) {
var r = new TinyBigInt( 0 );
var old_r;
var x = new TinyBigInt( 1 );
var quot_rem = [ new TinyBigInt( 0 ) ];
var a = n;
var b = k.mod( n );
if( b.is_negative ) b = b.add( n ); // (*1)
while( b.compare( 0 ) !== 0 )
{
old_r = r;
r = x.sub( quot_rem[ 0 ].mul( r ) ); // r := x - qr
x = old_r; // x := (old) r
quot_rem = a.div( b, true ); // a = quot_rem[ 0 ] * b + quot_rem[ 1 ] = qb + c
a = b; // a := b
b = quot_rem[ 1 ]; // b := c
}
if( a.compare( 1 ) === 0 )
return ( r.is_negative ) ? r.add( n ) : r ;
Append( "*** 1/k (mod n) doesn't exist:\n"
+ " k = " + k + "\n"
+ " n = " + n + "\n"
+ " gcd(k,n) = " + a );
return null;
}逆数があれば、それを返す。最大公約数が 1 でないときは逆数が存在しないので、メッセージを出して null を返す。(*1) は、入力された k が負でもいいようにする ための細工。
JavaScript の剰余演算子は、負の数については負の剰余を返し、最小非負剰余を返さない。ライブラリの mod() もこれと互換の動作になっている。
これを使って、点 P をそれ自身に足す関数を作る（mod n）。「点」の実装は、普通にやれば「プロパティー .x と .y を持つオブジェクト」だが、手っ取り早く長さ2の配列を使い、配列の1個目・2個目の要素がそれぞれ x-座標、y-座標だと約束しよう。点 は、長さ0の配列 [] で表すことにする。引数 a は、楕円曲線の式の1次の係数。
function Ell_double( P, a, n ) {
if( P.length === 0 ) return []; // P = O
if( P[ 1 ].compare( 0 ) === 0 ) return []; // P = -P
var x = P[ 0 ], y = P[ 1 ]; // x = x1 = x2, y = y1 = y2
var inv = Inv_mod( y.mul( 2 ), n ); // 1 / (2y)
if( inv === null ) return null;
var lambda = ( x.sqr().mul( 3 ).add( a ) ) // λ = (3x^2+a) / (2y)
.mul( inv ).mod( n );
var x3 = lambda.sqr().sub( x ).sub( x ); // x3 = λ^2 - x - x
var y3 = lambda.mul( x.sub( x3 ) ).sub( y ); // y3 = λ(x - x3) - y
return [ _mod( x3, n ), _mod( y3, n ) ];
}
function _mod( k, n ) {
var ret = k.mod( n );
return ( ret.is_negative ) ? ret.add( n ) : ret ;
}逆数が null になって λ が計算できない場合は null を返し、そうでなければ点の和を返す（§3 参照）。結果を返す前に、点の座標値については n で割った余り（最小非負剰余）に置き換えて [0, n − 1] の範囲に収まるように統一しておく（関数 _mod() は、この処理を行う）。
特別な場合として、P.length === 0 は P が [] であること、つまり P = を表す。一方、P[1] = 0 つまり y1 = 0 の場合、マイナスの点についての規約から P = −P（なぜなら y1 = −y1）。いずれの場合も、P + P = となる。
次に、点 P と点 Q を足す関数を作ろう。足し算の相手が の場合（*6）、および、ある点と「マイナスその点」を足す場合（*8）は、特別に処理する。P と Q が同一の点の場合（*7）は、上記の Ell_double を呼び出す。それ以外の場合、和の定義に従って、点 P + Q の座標を返す。ただし計算できないときは null を返す。a は、楕円曲線の式の1次の係数。
function Ell_add( P, Q, a, n ) {
if( P.length === 0 ) return Q; // P = O (*6)
if( Q.length === 0 ) return P; // Q = O (*6)
var x1 = P[ 0 ], y1 = P[ 1 ];
var x2 = Q[ 0 ], y2 = Q[ 1 ];
if( x1.compare( x2 ) === 0 ) { // if x1 = x2 (*7)
if( y1.compare( y2 ) === 0 ) // if y1 = y2, Q = P (*7)
return Ell_double( P, a, n );
else return []; // otherwise Q = -P (*8)
}
var inv = Inv_mod( x1.sub( x2 ), n ); // 1/(x1 - x2)
if( inv === null ) return null;
var lambda = y1.sub( y2 ).mul( inv ).mod( n ); // λ = (y1 - y2)/(x1 - x2)
var x3 = lambda.sqr().sub( x1 ).sub( x2 ); // x3 = λ^2 - x1 - x2
var y3 = lambda.mul( x1.sub( x3 ) ).sub( y1 ); // y3 = λ(x1 - x3) - y1
return [ _mod( x3, n ), _mod( y3, n ) ];
}2016年8月28日追記: (*7) では、座標値が「整数として等しいか」を判定している。理論上、この方法では、座標値が「整数として等しくないが、考えている法では合同」というケースが正しく処理されない。この記事のコードでは、この関数に渡される点の座標値を最小非負剰余に統一することにより、座標値が「合同なら整数として等しい」ことを保証して、この問題を回避している。
以上で「点の2倍」と「点の和」を計算する関数が準備された。これらを組み合わせると、「繰り返し2乗法」風に「点の m 倍」を計算できる。ただし途中のどこかで和が定義されず計算不能になった場合、戻り値 ret は null になる。
function Ell_mul( P, m, a, n ) {
var bin = m.toString( 2 );
var len = bin.length;
var ret = P;
for( var i = 1; i < len; ++i )
{
ret = Ell_double( ret, a, n );
if( ret === null ) break;
if( bin.charAt( i ) === "1" ) {
ret = Ell_add( ret, P, a, n );
if( ret === null ) break;
}
}
return ret;
} 準備ができたので、試しに B1 = 11000 として n = 2149 − 1 の因数を探してみる。Get_m は定義に基づき最小公倍数を返す関数で、これは別途実装しておく（ソースコード参照）。ここでは y2 ≡ x3 + ax + 1 (mod n) の形の楕円曲線を使い、曲線上の点として J = (0, 1) を選択する。
function ECM_Test( n, B1 )
{
var m = Get_m( B1 );
var J = [ new TinyBigInt( 0 ), new TinyBigInt( 1 ) ];
for( var a = 1; a < 100; ++a ) {
Append( "ECM: y^2 ≡ x^3 + " + a + "x + 1 (mod n), B1=" + B1 );
var K = Ell_mul( J, m, a, n );
if( K === null ) break;
}
}
var n = TinyBigInt.pow( 2, 149 ).sub( 1 );
ECM_Test( n, 11000 );「p − 1」法では手が届かなかった「ハロロ・コロニ・ロハコロロニハニ・イハサイコイ」がついに検出された！
和の計算の途中で「逆数がありません」というエラーを引き起こしたこの数こそ、まさに求めたかった n の因数だ。
上記のテストでは、11種類の楕円曲線について K = [m]J (mod n) が計算された。その一つにおいて逆数が定義されないケースが発生し、原因の数が報告され、K が null となってループが終了した。計算で使われた各曲線について、p = 86656268566282183151 を法とした場合の点の数（位数）は次の通り（これらの数は PARI/GP（）で計算された）:
いずれの曲線でも位数は大ざっぱに p だが、正確な位数はかなりランダムだ。結果として、位数を素因数分解したときの最大因数の大小も、曲線によって大きく異なる。p − 1 における「38兆」よりさらに大きいケースもあるし、逆にずっと小さいケースもある。特に a = 11 とした楕円曲線 y2 ≡ x3 + 11x + 1 (mod p) の位数は、10753 以下の素因数しか含まず、B1 = 11000 のとき、因数分解成功の条件が満たされる。
以上はあくまで原理を紹介するための実演であり、実用上、このままでは速度的な不満がある。テスト環境では曲線1個で5秒かかり、合計約55秒。同じく20桁の因数を持つ 2137 − 1 も分解できるが、その場合 a = 1016 まで続ける必要があり、ざっと1時間半もかかってしまう。
しかし成功に必要な B1 を大幅に減らせることは確かなので、効率を改善できれば実用性が高まる。
具体的に、標準形の楕円曲線の代わりにモンゴメリー形式を使うことで、計算が速くなる。さらにステージ2を併用することで、実効性が高まる。加えて、この例のようにメルセンヌ数が法になる場合、単純に mod を呼び出さずに剰余演算の実装を工夫することで、桁違いの高速化が実現される。
米国の数学者ピーター・モンゴメリー（ Peter Lawrence Montgomery）は「計算の鬼」だ。「これはすごい」と思わせるアルゴリズムをいくつも発見している。以下で紹介する計算法では、楕円曲線の標準形の代わりに、
という表示（モンゴメリー形式）が使われる。さらに、点の位置が「x-座標と z-座標を使い、y-座標を使わない」という独特な方法で表現される。
通常の座標（アフィン座標）の (x/z, y/z) と同じ意味で、x : y : z という比を考える。座標であることを明確にするため、この比を (x : y : z) と書いて、3個の要素をそれぞれ x-座標、y-座標、z-座標と呼ぶことにする。これは「比」なので、3個の座標値を定数倍（0倍以外）したり、0以外の定数で割ったりしても、表す点は変わらない。このような点の表し方は、同次座標と呼ばれる。斉次座標または射影座標と呼ばれることも多い。
z = 1 なら、これは普通の座標の (x, y) と同じだし、座標値を z で割り算すれば、いつでもこの形に戻すことができる（z の逆数が存在する限りにおいて）。
楕円曲線法の文脈において、このような座標系を使う主目的は処理の効率化。逆数の計算は面倒なので、割り算を先延ばしして比の形で座標を保持し、どうしても必要になった場合だけ割り算する。
楕円曲線の標準形を使う §4 の方法では、毎回の λ の計算において「分母になる数の、mod n での逆数が存在するか？」が問題の核心だった。モンゴメリーの方法では、同じ内容のチェックが曲線1個につき1回で済む。K = [m]J を計算した後で K の z-座標と n の最小公倍数を調べ、それが 1 でなければ、通常は分解成功となる。
同次座標は、モンゴメリー形式に限らず、標準形の楕円曲線でも利用可能。座標系をきちんと定義した場合、無限遠点 の定義も まともになる。すなわち、 = (0 : 1 : 0)。この記事では座標系をきちんと定義しないが、とりあえず「無限遠点は z ≡ 0 の点」と考えることができる。
モンゴメリーの計算法の大胆さは「y-座標を使わない」形になっているところ。「y-座標は必要ないので、計算しなくていい」と割り切っている。
この方法では、実質的に点の x-座標 x/z しか分からず、「点 P」と呼ばれるものは、実際には「点 P の x-座標の分子と分母」にすぎない。計算法の詳細は [Mon87] の 260–261 ページに記されている。P + Q の計算で P − Q を利用するところが特徴的。
メルセンヌ・ウィキの Formulas for addition and duplication にも説明がある。 Montgomery arithmetic も参考になるが、そこで 3M+2S
と記されているのは Zm−n ≡ 1 というショートカットを利用した場合の話で（補足3）、一般には 4M+2S となる。
2016年9月4日追記: モンゴメリーの計算法の公式の導出については、[Sut15] §11.9 で詳しく説明されている。
以下、JavaScript による実装例を示す。y-座標不要のため、点 (x : y : z) を長さ2の配列 [x,z] として実装している。コメントでは、これを (x::z) と書いている。2重コロンは単に「3項の比の真ん中が抜けたもの」で、スコープ解決演算子ではない。
function Mon_double_Me( P, A24, Me, e ) { // P = (x::z)
var Q = new Array( 2 ), u, v, s, t; // P+P = Q = (x'::z')
u = _modMe( P[ 0 ].add( P[ 1 ] ).sqr(), Me, e ); // u = (x+z)^2
v = _modMe( P[ 0 ].sub( P[ 1 ] ).sqr(), Me, e ); // v = (x-z)^2
Q[ 0 ] = _modMe( u.mul( v ), Me, e ); // x' = u * v
s = u.sub( v ); // s (=4xz) = u - v
t = _modMe( v.add( A24.mul( s ) ), Me, e ); // t = v + ((A+2)/4)s
Q[ 1 ] = _modMe( s.mul( t ), Me, e ); // z' = s * t
return Q;
}A24 は（モンゴメリー形式の楕円曲線の）2次の項の係数 A に関係する値で、(A + 2)/4 を表す。この値は、呼び出し側で計算される。
_modMe(x, Me, e) は、メルセンヌ数 Me = 2e − 1 について x mod Me を簡約する。x.mod(Me) より高速に動作し、0 以上 Me 未満の値を返す。実装は、前回の _modMp と同様。ここでは法がメルセンヌ数なのでメルセンヌ数で簡約しているが、これは楕円曲線法とは関係ない部分であり、一般には n を法とする簡約（還元）を行えばいい。
function Mon_add_Me( P, Q, P_minus_Q, Me, e ) { // P = (x1::z1), Q = (x2::z2)
var R = new Array( 2 ), u, v, s, t; // P+Q = R = (x3::z3), P-Q = (x0::z0)
u = _modMe( P[ 0 ].sub( P[ 1 ] )
.mul( Q[ 0 ].add( Q[ 1 ] ) ), Me, e ); // u = (x1-z1)(x2+z2)
v = _modMe( P[ 0 ].add( P[ 1 ] )
.mul( Q[ 0 ].sub( Q[ 1 ] ) ), Me, e ); // v = (x1+z1)(x2-z2)
s = _modMe( u.add( v ).sqr(), Me, e ); // s = (u+v)^2
R[ 0 ] = _modMe( P_minus_Q[ 1 ].mul( s ), Me, e ); // x3 = z0 * s
t = _modMe( u.sub( v ).sqr(), Me, e ); // t = (u-v)^2
R[ 1 ] = _modMe( P_minus_Q[ 0 ].mul( t ), Me, e ); // z3 = x0 * t
return R;
}P_minus_Q は、楕円曲線上の点 P − Q すなわち P + (−Q) を表す。点 P と点 Q を足したいときに、それらの「差」という第3の点が必要になるのは不便だが、ステージ1では、このことは大きな問題にはならない。下記の「はしご」を使うと、2点の差を一定に保ったまま計算を進められるからだ。
「2倍」と「和」のいずれの実装でも、点が であるケースをきちんと扱っていない。「和」の実装では P = ±Q のケースも無視されている。しかしこれらの「特殊事態」を厳密に扱うことは難しい（その気になれば P = ±Q になっていないかチェックできるが、速度的に重荷となる）。特殊事態が発生する可能性は低い上、万一発生して計算がおかしくなっても「その曲線では因数分解できない」という以上の害はないので、特殊事態については考えないことにしよう。
上記の「2倍」と「和」を組み合わせて「点の定数倍」を求めるアルゴリズムはモンゴメリーのはしごと呼ばれる。
function Mon_mul_Me( P, m, A24, Me, e ) {
var kP, k1P, Q; // kP = [k]P, k1P = [k+1]P, Q = [2k+1]P
var bin = m.toString( 2 );
var i, len = bin.length;
kP = P; // [k']P = [1]P (k'=1)
k1P = Mon_double_Me( kP, A24, Me, e ); // [k'+1]P = [k']P + [k']P = [2]P
for( i = 1; i < len; ++i )
{
Q = Mon_add_Me( k1P, kP, P, Me, e ); // [2k+1]P = [k+1]P + [k]P (k=k')
if( bin.charAt( i ) === "0" ) {
kP = Mon_double_Me( kP, A24, Me, e ); // [k']P = [k]P + [k]P = [2k]P (k'=2k)
k1P = Q; // [k'+1]P = [2k+1]P
} else {
kP = Q; // [k']P = [2k+1]P (k'=2k+1)
k1P = Mon_double_Me( k1P, A24, Me, e ); // [k'+1]P = [k+1]P + [k+1]P = [2k+2]P
}
}
return kP; // [k']P = [m]P (k'=m)
}for ループの間中、2点 [k]P と [k+1]P は次々と変化するが、それらの差は常に P に等しい。m を最初に計算して曲線ごとに Mon_mul_Me を1回だけ呼び出す場合、この「2点の差」は、掛け算 [m]J の対象となる最初の点 J に他ならない。
「モンゴメリーのはしご」はシンプルで美しいが、最速のアルゴリズムではない。高速化を追求した実装では、代わりに「リュカの鎖」と呼ばれるアプローチが使われる。
σ（シグマ）と呼ばれる整数（6以上）を経由して、計算に使う最初の点 J の座標 (x0, z0) や曲線の2次の項の係数 A に対応する値が決定される。
この方法は、日本の研究者スヤマ・ヒロミ（*）のアイデアに基づく。オーストラリアのリチャード・ブレント（Richard P. Brent（））は、[Bre86] において「位数が12の倍数になるように、スヤマが示唆した実装を行った」と述べている。スヤマの媒介変数表示とも呼ばれる公式は、[Bre99] の434ページ、[Zim06] の527ページに記されている。
「位数が12の倍数」ということは、一般の位数と比べると検出したい因数が1桁小さくなる ようなもので、一見地味だが実効性が高い（§7 参照）。1985年に [Bre86] の初版（**）が公開されたとき、スヤマはブレントに私信を送り、このような観察を伝えたらしい。
function Get_A2_and_J( sigma, J, Me, e ) {
var u, v, s, t, inv;
u = new TinyBigInt( sigma * sigma - 5 ); // u = σ^2-5
v = new TinyBigInt( 4 * sigma ); // v = 4σ
J[ 0 ] = _modMe( u.pow( 3 ), Me, e ); // x0 = u^3
J[ 1 ] = _modMe( v.pow( 3 ), Me, e ); // z0 = v^3
s = _modMe( v.mul( J[ 0 ] ).mul( 4 ), Me, e ); // 4u^3 * v = v * x0 * 4
inv = Inv_mod( s, Me ); // inv = 1/(4u^3 * v)
if( inv === null ) return null;
t = _modMe( v.sub( u ).pow( 3 ).mul( inv ), Me, e ); // t = (v-u)^3 / (4u^3 * v)
return _modMe( t.mul( u.mul( 3 ).add( v ) ), Me, e ); // A+2 = t * (3u+v)
}この関数は、配列 J の要素（つまり点 J の座標値）を書き込み、A+2 を返す。
function Mon_Test( Me, e, B1 ) {
var sigma, A2, J = new Array( 2 ), inv, A24, m = Get_m( B1 ), K, g;
inv = Inv_mod( new TinyBigInt( 4 ), Me ); // 1/4
for( sigma = 6; sigma < 400; ++sigma )
{
Append( "Mon_Test: sigma=" + sigma + ", B1=" + B1 );
A2 = Get_A2_and_J( sigma, J, Me, e ); // A2 = (A+2)
if( A2 === null ) break;
A24 = _modMe( A2.mul( inv ), Me, e ); // (A+2)/4
K = Mon_mul_Me( J, m, A24, Me, e ); // K = [m]J
g = _gcd( K[ 1 ], Me ); // (*2)
if( g.compare( 1 ) !== 0 ) {
Append( Me + " = " + g + " × " + Me.div( g ) );
break;
}
}
}
var e = 149;
var Me = TinyBigInt.pow( 2, e ).sub( 1 );
Mon_Test( Me, e, 11000 );この方法だと、曲線1個0.2秒で済む（テスト環境での実測値＝以下同じ）。1個5秒かかっていた最初のテストとは大違いだ。しかしラッキー・シグマが341と大きいので、曲線を300個以上試す必要があり、このままでは因数分解成功までに約75秒かかる。一方、M137 の20桁因数については、最初の方法では1時間以上必要だったが、Mon_Test では57秒で検出された。
ステージ2を併用すれば、どちらも10秒以内に分解される。
(*2) の行について、メルセンヌ・ウィキでは「x-座標との最大公約数を取る」としているが、z-座標との最大公約数を取る方が、ほんの少し因数の検出力が向上する。普通はどちらでも同じだが、「B1 の大きさがギリギリ足りている」場合、違いが生じる。
注: 初版では inv の計算を毎回ループ内で行っていたが、この計算をループに入る前に1回だけ行うように改めた（2016年9月4日）。これによって、理論上0.1%程度の高速化が見込まれる。ただし、実測では有意な速度差はなかった。「4で割る」処理には近道も考えられるが、これはテスト用のスケルトンなので、単純にそのまま計算している。
Get_m に関連する問題ついては、補足3参照。
ある σ について、ステージ1で因数が検出されなくても、続けてステージ2を行うことができる。多くの場合、成功はステージ1ではなくステージ2で起きるので、ステージ2の役割は大きい。
ステージ2では、ステージ1とほぼ同じ計算時間を使って、B1 のざっと100倍の大きさの数 B2 までを処理することができる。ただし、次の条件が満たされなければ ならない。
考えている楕円曲線 E について、mod p での点の数（位数）を #Ep として、それを素因数分解したとき:
…例えば B1 = 11000, B2 = 99000 とした場合、位数が
ならステージ2は成功するが、位数が
の場合、最大因数は第1の例より小さいにも かかわらず、ステージ2は成功しない。B1 と B2 の間に2個の因数があるからだ。
この制約はあるものの、素因数分解の結果は「特に大きい因数が1個。それ以外の因数はどれも比較的小さい」というパターンを持つことが多いので、ステージ2は実用上うまくいきやすい。
ステージ1においては、m が #Ep の倍数なら
(mod p)となって K = [m]J (mod n) が計算不能となり、分解成功するのだった。そうすると、ステージ2の成功条件が満たされる場合というのは、「ステージ1はほとんど成功だったが、m を #Ep の倍数にするには1個だけ因数が足りなかった。あと1歩、及ばなかった」という状況であることが分かる。足りなかった1個は #Ep の最大素因数であり、B1 と B2 の間に存在するある素数だ。それを α（アルファ）としよう。この α さえ m の計算に加わっていたら、[m]J = (mod p) となっていた。すなわち、今からでも遅くないから、次の計算をすればいい:
そうすれば、内部的に
(mod p)となって「楕円さま、ご乱心」が発生し、アフィン座標（普通の座標）なら「逆数がありません」エラーが起きて p が発覚するし、同次座標（射影座標）なら z-座標の因数として p があぶり出される。
以上がステージ2の動作原理であり、原理は単純だ。
素朴な方法でこれを実行するなら、B1 < si ≤ B2 を満たす全ての素数 si について [si]K (mod n) の計算を試みればいい。si のどれかが α なのだから、途中のどこかで成功条件が満たされる。
しかしナイーブなやり方では、速度的なメリットがない。[si]K を一つ一つ計算することは、掛け算の回数と手間という点で、次の計算とほぼ同じことだ:
このように計算する場合、係数をどのくらい まとめて掛けるかは別問題として、実質的に m そのものを再定義していることになり、単なる「ステージ1の規模拡大」になってしまう（B1 と B2 の間に2個以上の因数があっても構わないことになる）。
いくぶんまともなアプローチは、[si]K から次々と [si+1]K を計算することだろう。素数の間隔についての考察によれば、si と si+1 の差は一定範囲の正の偶数 d であり、前もってその範囲の全部の d について [d]K を計算してキャッシュしておけば、「ある点から次の点を導出する処理」は、その点と [d]K の足し算となる。
悪いアイデアではないが、以下では別の道筋を考えてみたい。
si が大きいとき、[si]K を一つずつ計算するのはコストがかかる。山の斜面にある家の一つ一つに向けて
具体的に、各素数 si を2個の整数の和として表すことを考える。
そのうち第1の整数は、あらかじめ決めてある大きめの定数 w（例えば w = 210）の倍数であり、si に近い値になるようにする。ここでは、これを vw とする。第2の整数 u は、si と vw の差に当たる部分で、小さな値を持つ。例えば、
と書くと、それを満たす v, u が si ごとに一意に定まる。ステージ2が成功する場合、仮定により si のどれか一つ（それを α と呼ぶ）について:
(mod p)ここで、左辺の [210v]K を点 G と名付ける。さらに、右辺の符号を変えた [u]K を点 H と名付ける。アフィン座標（普通の座標）で考えた場合、楕円曲線における「マイナスの点」の定義によれば、G と H は x-座標が等しい。つまり、点 G の x-座標 を Gx として、点 H の x-座標 を Hx としたとき、
となって、Gx − Hx は p の倍数。この関係を満たす si がどこかに あるのだから、各 si について上記と同様に点 G と点 H を定めて Gx − Hx を求め、それらを全部掛け合わせれば、結果（それを r とする）は p の倍数になる。ゆえに r と n の最大公約数を計算すれば p が得られる（r が q の倍数にもなってしまった場合は例外だが）。r を求める掛け算・その前提になる計算は、mod n でも同じ結果になる（この計算をする時点では p は不明なので、mod n を使うしかない）。
このうち、
の部分は、次のように計算される。まず「B1 を超える最小の素数」に対応する v を求めて、それを使って最初の G を計算する。その素数と同じくらいの大きさの素数に関しては、同じ G を使い回すことができる。それでは足りなくなったら（＝隣のバス停のエリアに入ったら）、G に [210]K を足して、結果をあらためて G とする。以下同様。点 G については、1ステップごとに [210]K 進む「巨人の足跡」だとイメージすることができる（1ステップごとの前進量が大きい）。
一方、
の部分は、数が限られているので最初に全部計算してキャッシュしておき、それをステージ2全体で使い回せばいいだろう。点 H の計算は [1]K 単位の小刻みなものであり、これについては「赤ちゃんの小さな足跡」とイメージすることができる。si は素数なので、この場合、H については u が奇数の場合だけ（もっと言えば、210 と互いに素な場合だけ）計算しておけば足りる。
ここで 210 という値は話を具体的にするための例であり、一般には何らかの定数 w ということになる。
点を同次座標（射影座標）で考えている場合、上記の「x-座標が等しい」という部分は、「x-座標と z-座標の比が等しい」となる。そうすると、上記の Gx ≡ Hx に当たる条件は:
つまり「巨人の点」と「赤ちゃんの点」の各組み合わせについて「x-座標と z-座標をたすきに掛けて引き算したもの」を計算して、それらの積を r (mod n) とすればいい。
ここまでの内容を整理する。一定の方法で定義される点の組 G, H のそれぞれについて Gx − Hx (mod n) を計算し、それらの積 r を求めると、r は探し求めている素因数 p の倍数になるのだった（ステージ2成功の前提が満たされた場合には）。ただし、これはアフィン座標（普通の座標）の場合であり、同次座標（射影座標）の場合、r の構成要素は「たすき掛けの差」になるのだった。以下ではアフィン座標で記述するが、「差」の部分を「たすき掛けの差」に置き換えれば、同次座標でもほとんど同じ議論となる。
さて、楕円曲線の定義から、点 H と点 −H は x-座標が等しい。ということは、Gx − Hx を計算したとき、それと等しい値として Gx + Hx も計算されている ことになる。
この観察は、重大な意味を持つ。
そもそも Gx − Hx (mod n) という値は、B1 と B2 の間の素数:
の一つ一つに対応して定まるものであり、一つの Gx − Hx を計算することは、「巨人ステップ」と「赤ちゃんステップ」を経由して
という点を計算し、その点が持つ意味を抽出することだ。ここで「意味」というのは、楕円曲線法の文脈における「意味」であり、はっきり言って、これらの点のほとんどは何の意味も持たず、1個や2個計算を飛ばしても通常は大勢に影響ないのだが、一つだけ「当たりの素数」 α に対応する「当たりの点」があり、その点を計算に加えることで因数分解が成立するのだった。「当たりの素数」がどれだか分からないので、しらみつぶしに計算している。
ところで上記の観察によれば、si に対応する点を処理した瞬間、
に対応する点も、自動的に計算に加わっている ことになる。すなわち si 自体が「当たりの素数」でなくても、たまたま si′ が「当たりの素数」に等しければ、結果は「当たり」となる。
si′ は素数 si の計算に付随して勝手に定まる整数であり、一般には素数ですらない のだから、こちら側で当たりを引く可能性は低そうだし、たまたま「裏で当たり」になって くれなくても、B1・B2 間の全素数を明示的に計算に含めている以上、当たりが出る楕円曲線では必ず「表で当たり」が出る。「裏のまぐれ当たりなど必要ない」とも思える。しかし、この現象を別の角度から考えてみよう。
次の二つの表現は、「赤ちゃんが歩く向き」が違うだけで本質的に同じ機能を持つ:
後者の表現の −105 ≤ u ≤ 104 について、u が負の場合には符号を反転させて正として計算しても、r は変わらない。u の符号を変えると si はもとの素数とは異なる値（一般には素数ですらない値）になるが、H = [u]K の x-座標は u の正負によらず一定であり、従って Gx − Hx は もとのまま だからだ。「u が負なら符号を変える」とすれば、「赤ちゃん」の歩数の範囲を半分に圧縮して、
において、全素数を処理できることになる。
このとき、210v ± u のプラスとマイナスがたまたま両方とも素数になっているケースにおいては、プラス側の素数とマイナス側の素数のどちらでも Gx − Hx は同じ値なのだから、210v + u についての1回の計算で素数が2個処理され、計算量が少し節約される。「1個の値段で2個買える」「一石二鳥」「桂で両取り」というような状況だ。
この現象を意図的に発生させられない だろうか？
例えば 6v ± 1 は、双子素数でない限り、両方素数にはならない。6v + 1 は素数だが 6v − 1 は素数ではない場合を考えてみる（例えば 127 & 125）。v を 1 減らして、この素数を 6v + 7 と表したらどうなるか。運が良ければ 6v − 7 は素数なので、この解釈によって「プラス・マイナスを両方素数」にすることができる（例えば 127 & 113）。それが駄目でも、6v ± 13 型のペア、6v ± 19 型のペアなど、いくつもの可能性が考えられる。
u の絶対値がある程度大きくなることを認めた上で、vw ± u のプラス・マイナスが両方素数になるように v, u の組み合わせを選ぶことを素数ペアリングと呼ぶ。うまく組み合わせれば、1度に2個ずつ素数が処理されるのだから、「個々の素数から意味を抽出する」という作業の効率が通常の2倍になり、2倍の高速化が見込まれる。
ステージ2には多くのバリエーションがあり、奥が深い。ここでは「赤ちゃんステップ・巨人ステップ」に基づく比較的シンプルな実装を行う。ただし、素数を1個ずつ処理するのではなく、素数ペアリングを併用する。
アイデアの大部分はモンゴメリーによる。ブレントは「誕生日パラドックス」方式という別のアプローチを考えている。
ペアリングの最も単純な例は「6v + 1 と 6v − 1 の少なくとも一方が素数のとき、6v + 1 に対応する G, H を計算する」。もう少し複雑な例は「30v ± u (u = 1, 7, 11, or 13) について 30v + u と 30v − u の少なくとも一方が素数のとき、30v + u に対応する G, H を計算をする」。いずれも積極的なペアリングではなく、たまたま ± の両方が素数になる場合、結果的に計算量が節約される。[Gaj06] には、そのような実装の疑似コードが記されている。出発点として、参考になるかもしれない。
積極的なペアリングの前提となるのは、ペアを表す式 vw ± u の u が w/2 より大きくなることを許容すること。そうすれば個々の素数について、それを表す v, u の組み合わせが複数存在するようになるので、選択の余地が生じる。目標は、範囲内の全素数をこの形で表し、なるべく ± がどちらも素数になるように（そして、なるべく同じ素数を2回以上含めないように）することだ。
具体的には、B1 < s ≤ B2 を満たす一つ一つの素数 s について、整数 v, u を用いて、次の設定でペアリングを行う:
ペアリングのために、モンゴメリーのアルゴリズム PAIR の必要部分を利用する。[Mon87] の251–3ページ参照。ただし umax は w の倍数だと仮定する（→ 補足4）。モンゴメリーは PAIR を「p − 1」法の文脈で紹介しているが、このアルゴリズムは楕円曲線法にも適用可能だ（同256ページ）。
モンゴメリーの待ち行列（キュー）たち
を文字通りに実装すると、例えば w = 210 の場合、Q[-209], Q[-208], …, Q[208], Q[209] となる。このような「負数のインデックス」は JavaScript のオブジェクトとして合法なので、以下では、この形をそのまま使う。本番の実装では、インデックスに w − 1 を足して普通の配列にした方が効率が良い（付属コードの PAIR_Test2 参照）。以上は「待ち行列をプロパティーとするオブジェクト」または「待ち行列を要素とする配列」であり、いずれの場合も、待ち行列自体は、単なる配列として実装される。
function PAIR_Test( B1, B2, w, u_max, v_min, v_max )
{
var Q = {}, q, PA, v, s, a, a1, u;
for( q = -w + 1; q < w; ++q ) Q[ q ] = []; // "queues"
PA = new Array( v_max + 1 ); // pairing_arrays (array of arrays)
for( v = v_min; v <= v_max; ++v ) PA[ v ] = [];
for( s = B1 + 1 + (B1 & 1); s <= B2; s += 2 ) // for each odd s in the range (B1,B2]
{
if( ! Is_odd_prime( s ) ) continue; // skip if s is not prime
a = Math.round( s / (2 * w) ); // a = round(s/2w) = floor(s/2w+0.5)
q = s - a * ( 2 * w ); // s = a(2w)+q, q < |w|
for(;;) {
if( Q[ -q ].length !== 0 ) { // IF Q[-q] is non-empty...
a1 = Q[ -q ].shift(); // take its front element a'
u = ( a - a1 ) * w + q; // and try this u
if( u <= u_max ) { // if u is not too big
PA[ a + a1 ].push( u ); // pair (a+a')w ± ((a-a')w+q) (*3)
break;
} else // u is too big
PA[ 2 * a1 ].push( Math.abs( q ) ); // force-pair (a'+a')w ± |q| (*4)
} else {
Q[ q ].push( a ); // ...ELSE push a into Q[q]
break;
}
}
}
for( q = -w + 1; q < w; ++q ) { // empty the queues
while( Q[ q ].length !== 0 ) { // this queue is still not empty
a1 = Q[ q ].shift(); // remove the front element a'
PA[ 2 * a1 ].push( Math.abs( q ) ); // force-pair (a'+a')w ± |q| (*5)
}
}
return PA;
}Is_odd_prime(k) は、奇数 k について、素数か素数でないかの真偽値を返す。この関数は別途実装しておく。
PA（すなわち pairing_arrays）は「ペアリング用の配列」の配列。
どの素数 s も、2w で割った結果 s = a(2w) + q に従って処理される。ただし q は [−w+1, w−1] の範囲にある整数で、負になることもある（絶対的最小剰余）。s が任意の整数なら q = −w になり得るが、s は素数なので、そのパターンは生じない。
ペアになる相手が決まるまで s は待機する。待機は、待ち行列 Q[q] の末尾に a を追加するという形で行われる（a と q の組み合わせで s が特定される）。
ただし、ペアになる相手が既に存在している場合には、待機は必要ない。すなわち…
もし s が出現した時点において、（インデックスの符号が反対の）待ち行列 Q[−q] に何らかのデータが存在している場合、先頭の要素を取り出して a′ とする。s′ = a′(2w) − q は以前から待機していた素数であり、このとき
は2素数 s, s′ を表す。その場合 v = a + a′, u = (a − a′)w + q と置くと vw ± u は素数ペアなので、条件 u ≤ umax が満たされれば、s と s′ は直ちにペアになることができる。条件が満たされる場合、ペア登録として u を配列 PA[v] に追加する。これが (*3) の意味だ。s′ は s より前の素数なので、s より小さい。このことから u ≥ 1 となる。
一方、条件が満たされない場合、s についてはまだ可能性があるものの、s′ とペアになり得る「新しい」素数は もはや現れない（これ以上待っても、u はますます大きくなるばかり）。その場合、s′ = a′(2w) − q と a′(2w) + q を強制的にペアにしよう（後者は一般には素数ではなく、素数だとすれば既に別の相手とペアになっているか、または範囲外）。このときには2個の整数:
がペアになるので、v = 2a′ として、u = | q | を配列 PA[v] に追加する。これが (*4) の意味だ。
強制ペアリングによって a′ に対応する素数 s′ を処理した場合、まだ Q[−q] が空でなければ、先頭の要素を取り出して あらためて a′ とし、それが表す素数をあらためて s′ として、上記と同様に s, s′ のペアリングを試みる。s とペアになる条件を満たす相手が現れないまま Q[−q] が空になってしまった場合、または最初から Q[−q] が空の場合には、s は待機しなければならない（待ち行列 Q[q] の末尾に a が追加される）。
…メイン・ループは以上の繰り返しになる。メイン・ループ終了後、全部の待ち行列が点検され、待機中の素数が残っていれば (*4) と同様に強制ペア処理される。これが (*5) だ。
配列 PA のインデックス v は、
の v であり、インデックス vmin 未満の要素は定義されない（若干メモリー的に無駄がある）。PA[v] は、その v に関して、ペアリング vw ± u で使われる正の整数 u の配列。
例えば B1=11000, B2=B1*60, w=210, u_max=8*w という設定において
PA[100] = [17, 19, 61, 101, 211, 467, 523, 611, 673, 751, 817, 839, 1247, 1291, 13]となるが（配列の中身はソートされていない）、その意味は「210v ± u の v が 100 のとき、u = 17, 19, 61, … を使ってペアを作ればいい」。具体的なペアは
となり、15組30個の数はほとんど全て素数になっている。ただし、最後の1組だけは
となって、強制ペア処理の結果の合成数 20987 = 31 × 677 を含む。
この例では、モンゴメリー形式の同次座標（射影座標）が使われている。w が奇数の2倍であることを前提とした実装。
function Stage2_Test( K, pairing_arrays, w, u_max, v_min, v_max, A24, Me, e )
{
var baby_steps, K2, u, S, G, G_minus_S;
var r, v, Gx, Gz, arr, j, H, GxHz, HxGz, old_G;
baby_steps = new Array( u_max + 1 );
baby_steps[ 1 ] = K; // [1]K
K2 = Mon_double_Me( K, A24, Me, e ); // [2]K
baby_steps[ 3 ] = Mon_add_Me( K2, K, K, Me, e ); // [3]K
for( u = 5; u <= u_max; u += 2 ) // [5]K, [7]K, [9]K, ... (*9)
baby_steps[ u ] = Mon_add_Me( baby_steps[ u - 2 ], K2, baby_steps[ u - 4 ], Me, e );
S = Mon_double_Me( baby_steps[ w/2 ], A24, Me, e ); // [w]K : one Step of G (the Giant)
G = Mon_mul_Me( S, v_min, A24, Me, e ); // G := [v_min]S = [v_min*w]K
G_minus_S = Mon_mul_Me( S, v_min - 1, A24, Me, e ); // G - S = [v_min]S - S = [v_min-1]S
r = new TinyBigInt( 1 ); // r := 1
for( v = v_min; v <= v_max; ++v )
{
Gx = G[ 0 ], Gz = G[ 1 ]; // coordinates of G
arr = pairing_arrays[ v ]; // array PA[v]
for( j = 0; j < arr.length; ++j ) {
u = arr[ j ]; // how many steps does the Baby take
H = baby_steps[ u ]; // H := [u]K
GxHz = _modMe( Gx.mul( H[ 1 ] ), Me, e ); // Gx * Hz
HxGz = _modMe( H[ 0 ].mul( Gz ), Me, e ); // Hx * Gz
r = _modMe( r.mul( GxHz.sub( HxGz ) ), Me, e ); // r := r(Gx * Hz - Hx * Gz)
}
old_G = G;
G = Mon_add_Me( G, S, G_minus_S, Me, e ); // G := old G + one Step
G_minus_S = old_G;
}
return r;
}前半では「赤ちゃんステップ」テーブルが作成される。赤ちゃんの1回の歩数は K を単位に 1 以上、umax 以下なので、その範囲の u について、赤ちゃんの座標 [u]K を計算してキャッシュしておく。ただし、u と w は互いに素であり、ここでは w が偶数だと仮定しているので、奇数の u についてのみデータを作れば十分。
さらに「巨人の一歩」に当たる S = [w]K が計算され、それを使って巨人の初期位置 G および一歩手前の位置 G − S が決定される。
後半の処理は §6 参照。巨人の位置 G の一つ一つについて、（その位置に対応する）ペアリング用の配列 arr から配列の長さ回、赤ちゃんの歩数が読み出され、赤ちゃんの位置 H は、歩数に対応する計算済みの座標となる。G と H の組み合わせの一つ一つは、B1・B2 間のどれかの素数（大抵はペアリングにより2個の素数）に対応している。配列 arr の全要素が処理されると、巨人が1歩前進し、G と G − S が更新される。
2016年9月11日追記: 本来 w と互いに素な u についてのみ [u]K があれば十分だが、モンゴメリー形式では P + Q の計算に P − Q が必要。例えば、[13]K, [17]K, [19]K だけを生成したいとしても、[17]K + [2]K = [19]K を計算するとき [17]K − [2]K = [15]K の座標が必要になる。この問題を手っ取り早く回避するため、(*9) では奇数の u について一律に [u]K を生成している。この部分については、工夫すればもう少し計算量を節約できる。
§5 の Mon_Test にステージ2を追加する。
function Mon_Test2( Me, e, B1 ) {
var sigma, A2, J = new Array( 2 ), inv, A24, m = Get_m( B1 ), K, g, r;
var B2 = 60 * B1;
var w = 210;
var u_max = 8 * w;
var v_min = Math.floor( B1 / w );
var v_max = Math.ceil( B2 / w );
var pairing_arrays = PAIR_Test( B1, B2, w, u_max, v_min, v_max );
inv = Inv_mod( new TinyBigInt( 4 ), Me ); // 1/4
for( sigma = 6; sigma < 400; ++sigma )
{
Append( "Mon_Test2: sigma=" + sigma + ", B1=" + B1 + ", B2=" + B2 );
A2 = Get_A2_and_J( sigma, J, Me, e ); // A2 = (A+2)
if( A2 === null ) break;
A24 = _modMe( A2.mul( inv ), Me, e ); // (A+2)/4
K = Mon_mul_Me( J, m, A24, Me, e );
g = _gcd( K[ 1 ], Me );
if( g.compare( 1 ) === 0 ) {
r = Stage2_Test( K, pairing_arrays, w, u_max, v_min, v_max, A24, Me, e );
g = _gcd( r, Me );
}
if( g.compare( 1 ) !== 0 ) {
Append( Me + " = " + g + " × " + Me.div( g ) );
break;
}
}
}for ループの前で、全部の曲線に共通の設定・初期化を行う。ここでは B2 を B1 の60倍として、umax = 8w = 1680 としている。さらに、ループ内において「ステージ1が駄目なら続けてステージ2を行う」という処理が追加された。
M149 について Mon_Test2( Me, e, 11000 ) を実行すると、今度は8秒で分解成功した。ラッキー・シグマは25だった。シグマ25の曲線の（mod p における）群の位数は、次の形になる:
これは「ステージ2の成功条件」で取り上げた例そのもので、B1=3200, B2=83000 程度で足りることが分かる。
M137 について、同じ Mon_Test2( Me, e, 11000 ) を使うと11秒で分解された。ラッキー・シグマは33で、群の位数は:
これはステージ2の効果を示す好例だろう。ステージ1だけだと B1 をざっと50万にする必要があるが、ステージ2を使った場合、B1=10000 とその50倍程度の B2 で軽快に対応できる。
B1=50000 を使い Mon_Test2( Me, e, 50000 ) とすると、ラッキー・シグマが8になり、M137 は5秒半で分解される。群の位数は:
使われた設定 B1=50000, B2=60*B1 は成功条件を満たしている。
…いずれも、分解が速いのはステージ2のおかげだ。ステージ2の成功条件は比較的成り立ちやすいので、少ない個数の曲線を試すだけで済む。さらに「赤ちゃん・巨人」と素数ペアリングが、ステージ2そのものの高速化に寄与している。上記のコードでは行っていないが、大きめの B1 に対しては、ペアリング設定も大きめにする方がいい（例えば w = 2310, u_max = w*4）。
どの位数も必ず 22 × 3 を約数として含んでいる。これがスヤマ・チューニングの効果で、「お金を崩す」ように位数を「崩す」働きを持つ。位数の「総額」が同じでも大きい素因数が発生しにくくなり、小さめの B2 でもアルゴリズムが成功しやすい。
B2 が大きければ大きいほど、より大きな「位数の素因数」に対応できる。しかし B2 が大き過ぎると曲線1個当たりの処理が遅くなり、トータルでの効率は悪化する。文献上、B2 の相場は B1 の50～100倍程度。目安として、ステージ1にかける時間とステージ2にかける時間がだいたい同じになるのが妥当とされる。ステージ2で強力な高速化を行っている場合、100より高い倍率が妥当になることもある。
メルセンヌ合成数を対象に上記とほぼ同じ実装を試したところ、20～30桁くらいの因数は、数分から十数分で検出されることが多いようだ（付録B 参照）。所要時間を気にしなければ、30桁くらいまでの因数は検出される。数時間～数日を費やすなら、30～40桁の因数も検出される可能性がある。JavaScript 上の簡易的な実装にしては、そう悪くないだろう。
この方法を使っても、現状、100桁や200桁の因数は検出できない。桁数が増えると時間がかかり過ぎて実用にならない…という性質は、素朴な試行除算や「p − 1」法と変わらない。「任意の整数を因数分解する速い方法があるのか」という根本的問題は、依然として未解決だ。
楕円曲線法（ECM: Elliptic Curve Method）は「p − 1」法の発展形だが、常に「p − 1」法より優秀とは限らず、状況によっては「p − 1」法の方が高速であることも珍しくない。
それでも ECM は、検出できる因数の範囲をある程度、広げてくれる。「p − 1」法では20桁くらいの因数が限界かもしれないが、ECM なら30桁くらいまでは いける。条件によっては40桁以上も可能だろう。「p − 1」法では位数がスムーズでなく（＝小さな因数だけの積ではなく）分解が難しい場合でも、ECM なら、位数がスムーズな群に出会うまで、曲線を変えて何度でも再試行できる。
この記事を書いた直接のきっかけは、前回の「[JS] メルセンヌ数の分類と分解」において、M149 が分解できなかったこと。メルセンヌ数を扱ったのは、もともとは任意精度ライブラリの動作テスト・実演のためだった。
この記事の前半と同等の内容（楕円曲線の標準形＆アフィン座標のステージ1）は、一般的な教科書にも記されている。モンゴメリー形式のステージ1は普通の教科書の範囲外だろうが、モンゴメリー自身の [Mon87] はオンライン公開されていて、自由に読むことができる。解説しているウェブページも複数ある。
多くの人が感じること だろうけど、「p − 1」法を知ってから取り組むと、楕円曲線法のステージ1は、あっけないほど簡単だ。ステージ1だけでも M149 は分解可能であり、大きな達成感が得られる。分解できないと、この20桁の因数が「38兆を要求する無法者」のようにも思えてくるが、分解できるとなると、全桁暗記するくらい愛着が湧く。
一方、楕円曲線法のステージ2については、ウェブ上で検索しても、見つかるのは素朴な実装と高度な話題の両極端で、中間が分からず、初めは手探り状態だった。素数の間隔を利用する実装なら「p − 1」法のステージ2の焼き直しで済むのだが、その先はどうなるのか…。ステージ2の原点に当たる [Bre86] [Mon87] だけでも内容は多岐にわたり、なかなかピンとこない。奥が深いから面白いのだが、それは難しいということ でもある。粘り強く、試行錯誤や実験を繰り返しながら、少しずつ理解できる範囲を広げていくしかない。[Gaj06] は標準版ステージ2を最小構成的に提示していて、最初の手掛かりとして大きな助けとなった。
この記事の話題の多くには、発展の余地がある。例えば「モンゴメリーのはしご」を「リュカの鎖」で置き換えること。「誕生日パラドックス」版ステージ2を試すこと。
素数ペアリングは「p − 1」法のステージ2にも適用可能だ。普通は「p − 1」法のテクニックを楕円曲線法に応用するのだが、素数ペアリングは逆順でやった方が自然だろう。この道順で進むと、ブレント・スヤマ拡張も自然と視野に入り、ステージ2でモンゴメリー形式から標準形に戻りたくなる理由も見えてくる。
「整数のようなもの」を座標値とする楕円曲線（きちんと言えば「n を法とする剰余類が作る有限環」上の楕円曲線）は、普通の意味での「曲線」ではない。「碁盤の上に散らばった碁石」とイメージしてもいいだろう。
碁石はランダムに散らばっている わけではなく、x-座標ごとに、
この違いは、「ある x を選んで固定したとき、曲線の定義式 y2 ≡ x3 + ax + b (mod n) を満たす y が存在するか しないか」という違いに対応する（右辺が平方剰余か非剰余か）。
この他、盤外の無限の彼方には という点がある。
n が素数の場合、任意の2個の碁石（同じ石を2回選んでもいい）に対して、それらの和と呼ばれる第3の碁石が矛盾なく定義される。n が素数の場合というのは、有限環上の楕円曲線の中でも特別なケースに当たり、有限体上の楕円曲線と呼ばれる。
楕円曲線法の厳密な扱いは難しいが、「曲線上の勝手な点を整数倍する」という本質は単純明快だ。ところが楕円曲線法は、必要以上に（というより中途半端に）小難しい表現で説明されることがある。
楕円曲線法の説明で往々にして不透明になっている部分について、（厳密に定式化するわけでは ないが）一歩踏み込んでみたい。
p と q が異なる素数で、n = pq とする。
楕円曲線法が成功する場合、[m]J = (mod p) かつ [m]J ≠ (mod q) となるが、そのとき [m]J (mod n) が計算不能になる（定義されない）。
「p − 1」法が成功する場合、同様に Jm ≡ 1 (mod p) かつ Jm ≢ 1 (mod q) となるが、Jm (mod n) が計算不能になる わけではない。
この違いは、どうして生じるのか？
§3 では楕円曲線 E: y2 ≡ x3 + 4x + 1 (mod 7) を調べ、この曲線上にある点の集合が、 を含めて
}であることを示した（5個の点が存在）。同じ曲線を mod 11 で考えると、曲線上には9個の点がある:
}同様に mod 77 での点の集合 E77 を考えることができる（33個の点が存在）。
今 E7 の点 A = (4, 2) と E11 の点 B = (5, 6) を選んだとすると、それらに対応して E77 の点 C = (60, 72) が定まる（60は7で割ると4余り、11で割ると5余る。72は7で割ると2余り、11で割ると6余る）。言い換えれば、C = A (mod 7) かつ C = B (mod 11)。一般に、E7 と E11 から一つずつ 以外の任意の点を選んだ場合、それらに対応する E77 の点が一意に定まる。逆に、E77 の任意の点（ を含む）に対して、E7 の点と E11 の点の組が一意に定まる。
「p − 1」法の場合、中国剰余定理により、群 (Z/nZ)* は群 (Z/pZ)* × (Z/qZ)* と同型であり、後者において成分ごとに計算できるものは、前者において mod n で計算でき、その逆も成り立つ。要するに、計算不能になることはない。
これと同様のことが、楕円曲線の点の群でも成り立つだろうか？
環 Z/pZ 上、および環 Z/qZ 上においてある一つの楕円曲線 E を考え、この曲線上の点（ を含む）の集合をそれぞれ E(Z/pZ) と E(Z/qZ) とする。E(Z/pZ) と E(Z/qZ) は、どちらも点の加法に関して群を成し、直積
も（成分ごとの）加法に関して群を成す。加法を持つ点の集合
を同様に定義する。もし仮に中国剰余定理の楕円曲線バージョンが成り立つなら、次の性質を持つ同型写像 f: G → S が存在して、S は 群になるはずだ:
ここで α, β は G の任意の2元。
上記の具体例では、G の位数は 5 × 9、S の位数は 33 で、前者の方が多いので1対1にならないことは明らかだが、状況を分析してみよう。
元 α = (Ap, Aq) の成分がどちらも でないなら、対応する S の点 f(α) = An が自然に定義される。単に An = Ap (mod p), An = Aq (mod q) になるように点 An の座標を選択すればいい（そのような選択が可能で一意的であることは、中国剰余定理によって保証されている）。Ap = Aq = の場合にも、f(α) = とすればいいだろう。しかし Ap, Aq の一方が で他方が でない場合、以下で述べる理由により f(α) は定義されない。言い換えると、f(α) と f(β) が定義されていても、f(α) + f(β) が定義される保証はない。すなわち、集合 S は群 G と同型でないばかりか、群でさえない。そのため「p − 1」法とは状況が異なり、条件を満たす写像 f は存在しない。
実際、楕円曲線の点の群では、和を定義するために、点の座標の四則演算が必要とされる。ところが環には乗法について可逆でない元が含まれており、点の和の計算でこの元が λ の分母になった場合、問題が生じ得る。
このうち λ の分母が ≡ 0 になりそうなケースについては、点の和の定義そのものによって、問題が回避される。すなわち「符号が反対の2個の点」を足し合わせる場合、定義により、和は になり、λ の計算は必要とされない。
λ の分母 μ が ≡ 0 でなく、かつ法と互いに素でない場合には、μ−1 が定義されず、従って λ が定義されない。法が素数なら そんなことは起きないが、mod n ではこの状況が生じ得る。それは gcd(μ, n) ≠ 1 の場合であり、言い換えれば μ ≡ 0 (mod p) かつ μ ≢ 0 (mod q) の場合、または μ ≢ 0 (mod p) かつ μ ≡ 0 (mod q) の場合だ。「G の成分の一方だけが 」と言い換えることができる。ただし、モンゴメリーの方法で割り算を先延ばししている場合には、実質的に計算不能に陥っていても、見掛け上は計算を続行できる。
…座標値レベルにおいて、 の意味は「μ ≡ 0 の乗法逆元は不要」ということ、 の意味は「μ ≢ 0 が乗法逆元を持たないと困る」ということ。法が素数なら は発生せず、環 Z/pZ の代わりに有限体 Fp を考えれば（同様に、Z/qZ の代わりに Fq を考えれば）話は丸く収まる。対照的に、環 Z/nZ は、どう頑張っても体にならない。これが「計算不能」の核心だろう。
の問題が起きたとき、確かに計算は破綻するものの、gcd(μ, n) が n の非自明な約数になっている のだから、楕円曲線法の目的は達成される。すなわち:
であるもの」については、対応する元が S にはない。その破れ目を踏んだ場合、S 側で和が定義されないが、そのときには探している因数が分かり、因数分解に成功する。本文中では Get_m を使って m を明示的に（一気に）計算して、点 J に対する掛け算は、1個の曲線につき1回で済ませている。つまり K = [m]J をそのまま計算している。
一般には m を明示的に計算せず、その構成要素の素数べき（または素数）を少しずつ点に掛け算する方がいい（後述の Tiny_ECM の「モード3」）。明示的に計算すると m は極めて大きな数になり、メモリー不足になったり動作が緩慢になったりする恐れがある。
ただし B1 が比較的小さい場合（例えば B1=50000）、m を明示的に計算した方が速いようだ。「少しずつ掛け算」方式では曲線ごとに暗黙に m を再計算することになるが、「一気に掛け算」なら最初に1回だけ m を計算しておいて、それを全部の曲線に対して使い回せる。これがプラスに作用するのだろう。ブレントもモンゴメリーも m を明示的に計算する必要はないと示唆しているが、当時と現在ではメモリー環境が全く違い、議論の前提が異なる。
モンゴメリー形式の曲線において「一気に掛け算」する場合、点 J = (x0::z0) を J′ = (x0/z0::1) と書き直しておけば、モンゴメリーのはしごの1段ごとの掛け算回数（mul または sqr の呼び出し）が11回から10回に削減される（Tiny_ECM の「モード2」）。実際には、σ が小さい場合、このショートカットは逆効果になる可能性が高い（掛け算回数は減るものの、J′ の x-座標は J の x-座標より桁違いに大きいため、かえって計算量が増えてしまう）。
B1 も sigma も小さい場合、「一気に掛け算するがショートカットはしない」のが良いようだ（Tiny_ECM の「モード1」）。
上記は楕円曲線法の場合であり、「p − 1」法の場合、m を明示的に計算するうまみは少ない（キャッシュしても使い回せないので）。
PAIR についてモンゴメリーは、umax ≥ w/2 で良いかのような書き方をしている（[Mon87] の251ページ）。それどころか umax = ⌊w/2⌋ が可能であるかのように扱っている（同253ページ）。
原理的にはそれは正しい主張だが、具体的なアルゴリズム PAIR に関する限り、実際には umax ≥ w − 1 でないと条件 u ≤ umax が満たされないと思われる。なぜなら、強制ペア処理のときには u = | q | だが、| q | の最大値は w − 1。
PAIR は umax = w/2 でも動作するが、その設定だと、要求仕様 u ≤ umax を満たさないペアが生成される可能性がある。その意味において、253ページの表2の1段目は、条件に基づく正確なカウントではない（PAIR が生成するペア数の正確なカウントではある）。
別のアルゴリズムを使えば、umax = w/2 という設定でも仕様を満たすペアリングは可能だが、その場合、全ての素数は一意的に vw ± u の形で書かれることになり、選択の余地がない。積極的なペアリングを行うためには umax がある程度大きいことが必要だ。
従って、umax = w/2 で PAIR が正常動作しないという事実は、実用上問題にならない（ペアリングによる高速化という観点からは、仮に正常動作するとしても、その設定を使う意味がない）。条件を umax ≥ w に変えれば問題は消滅する。
この記事の初版（2016年8月14日）では以下が使われた。
_26_.js: 任意精度整数演算ライブラリ、バージョン0.3（2016年8月7日版）。ecm.js: 本文で紹介されている全部のコードを含む（2016年8月14日版）。ecm.html: デモ用コードを呼び出す HTMLソースの例。呼び出しの大部分は、初期状態ではコメントアウトされている。試したい部分について1回に1個だけコメントアウトを外し、実行後またコメントアウトするといいだろう。テスト環境において1分以上かかるデモもある。JavaScript が遅い環境（ある種のブラウザなど）では、何十分もかかる場合があるので注意。本文（テキスト・画像）とソースは、全てパブリックドメイン。
2016年9月4日: ecm.js バージョン2を公開した。旧バージョンでは、Mon_Test, Mon_Test2 において定数 1/4 (mod Me) を毎回ループ内で計算していた。新バージョンでは、これをループに入る前に1回だけ計算するように改めた。それ以外に実質的な変更はない。本文中のコードも同様に更新（詳細）。
Tiny_ECM( Me, e, B1, B2, sigma_min, sigma_max, cofactor, mode ) は、付属ソースコード内の関数。本文で述べた方法を基に「メルセンヌ数に対するプチ楕円曲線法」を行う。第3引数以降は省略可能（省略された場合、デフォルト値が使われる）。
M193 を因数分解。
var e = 193;
var Me = TinyBigInt.pow( 2, e ).sub( 1 );
Tiny_ECM( Me, e );0.4秒で分解された。楕円曲線を持ち出すまでもない小さい素因数だが、これはただのテスト。問題は、この余因数の分解だ。
var e = 193;
var Me = TinyBigInt.pow( 2, e ).sub( 1 );
var cofactor = Me.div("13821503");
Tiny_ECM( Me, e, 50000, 50000*60, 6, 300, cofactor );ラッキー・シグマは121。検出された23桁の因数と余因数は、いずれも素数。所要時間は4分20秒だった（モード2では4分19秒)。B1=11000, B2=B1*60 という設定でも、sigma=415 まで粘れば分解に成功する（所要時間3分45秒）。
B1=50000, B2=B1*60 という設定を使って、M263 の第3因数（27桁）を検出。所要時間2分27秒。
M593 の4027
PDF]: 楕円曲線法の発見から2006年までの約20年間において、どのようなアルゴリズムの発展があったか概観している。 論文 / スライド]: シンプルなステージ2の疑似コードを含む。ステージ2入門となり得る。 PDF版 の方が読みやすい）。Integer Factorization Using Elliptic Curves から楕円曲線法の話になる。この記事の前半の内容について、数学的観点から分かりやすく記されている。明快な標準形ステージ1入門。
2016年10月2日追記: 上記の改訂版に当たる Elementary Number Theory: Primes, Congruences, and Secrets が公開されていて、誤植なども随時修正・更新されている。本文中で参照した逆数計算のアルゴリズムは、新バージョンでは §2.3.1 に当たる。
標準形のステージ1を扱った日本語の教科書もいくつかある。今回資料として使ったわけではなく、このような場所で間接的に宣伝するスタイルを好まないので、具体的な題名・出版社名などは書かないけれど、目次を見て「楕円曲線を使った因数分解」というようなセクションがあれば、たぶんステージ1の紹介だろう。その数ページのために本を買うのも無駄なので、図書館で探したらいい のかもしれない。
追記: ウェブ上で読める日本語の資料として、楕円曲線法に関しては [Izu99] があり、楕円曲線入門としては [Izu13] がある（下記「追加の参考文献」参照）。
情報が英語中心なのは善しあしで、帝国主義のようなものだが、母語至上主義だと [Len87] がオランダ語、[Gaj06] がポーランド語などとなって、ますます読むのが大変になる。
この記事自体、日本語が分からない人にとっては嫌がらせ のようなもので、客観的に考えると、この内容に興味を持つ人が日本語圏にたくさん いるとも思えず、メルセンヌ・ウィキに英語で書いた方が、より多くの人に役立って社会貢献になるかもしれない。とはいえ、多様性も大切だろう。
記事公開（2016年8月14日初版）の後に発見した資料:
[Izu99] 伊豆 哲也 (1999): 楕円曲線法の高速化について: ECM についての日本語の資料（ウェブ上で自由に読めるもの）としては、現在（2016年）これが最良かもしれない。素数ペアリングを別にすれば、この記事で取り上げた話題のほとんど全てが、簡潔に解説されている。ステージ2のバリエーションの説明を含む。伊豆先生は「伊豆・高木アルゴリズム」の発見者でもある（* / **）。これは、モンゴメリーの計算法（XZ座標）の一般化に当たる。
注: 基本の標準版ステージ2では、B1・B2 間の素数を抜き出して利用する。伊豆が紹介している「モンゴメリーのステージ2」は、それとは異なるバリエーション。この点は、[Tak95] も同様。
[Cra05] Richard Crandall; Carl Pomerance (2005): Prime Numbers: A Computational Perspective (Second Edition): コンピューター科学の専門家と著名な数論学者により、共同執筆された本。理論の簡略な説明と、多くの疑似コードから成る。モンゴメリーの計算法、楕円曲線法ステージ1・ステージ2の紹介を含む。ステージ2の実装では、素数の間隔を利用するシンプルなアルゴリズムをベースに、高速化が行われている。
[Izu13] 伊豆 哲也 (2013): 楕円曲線暗号入門（2013年度版）: 楕円曲線「暗号」と楕円曲線「法」は分野が違うが、入り口は共通。それが丁寧に解説されている。書庫内のファイル名は多分 Shift_JIS で、環境によっては文字化けするが、ファイルの中身は問題なく表示できた。
[Sut15] Andrew Sutherland (2015): Lecture Notes (MIT OpenCourseWare: Elliptic Curves, Spring 2015): このコースは、実装を視野に入れつつも「実用上動けばいい」というのではなく、数学的厳密性を重視しているようだ。これらの講義ノートには、他の場所にはあまり書かれていない ような貴重なヒントがいろいろ書かれている。最新（2010年代）の成果も紹介されている。最初に読む入門としては深過ぎるかもしれないが、文章は気さくで明快。講義ノート11（ PDF）のうち、§11.5 が「p − 1」法、§11.6 以下が楕円曲線法ステージ1。モンゴメリー曲線・モンゴメリーの計算法について詳しく説明されている。問題セット6（ PDF）には、ステージ2に関する理論が含まれている。
…この他、楕円曲線法に関する日本語の資料として、次のものが見つかった:
[Tak95] 高橋 大介; 鳥居 泰伸; 湯淺 太一 (1995): SIMD型超並列計算機における素因数分解: モンゴメリー形式を使ったステージ1とステージ2について、実装の細部にまで踏み込んで考察・実践している。ステージ2の「赤ちゃん・巨人」では、素数判定をせず全部掛け算する作戦を使っている。
[Izu00] 伊豆 哲也 (2000): 素因数分解に適した楕円曲線の生成法 [ PDF]: スヤマのアイデアとその発展について、明晰な文章で論じている。
[Koi01] 小池 慎一; 山住 富也 (2001): 楕円曲線法による素因数分解に関する実験的考察: 学生が書いた実験レポート。この記事と同じく、M137 などの分解を試みている。
[Mor15] 森下 拓也; 趙 晋輝 (2015): 擬似的2次拡大環上の楕円曲線法: 実装は示されていないが、スヤマ系の実験を行っている。エドワーズ曲線（比較的新しい話題）に言及している。
の両方をゼロと呼んでいる部分があったので改めた。Ell_add のテキストとコードの対応を番号で示し、(*7) に注釈を付けた。 (2) 表現の細部の調整（約15カ所）。 (3) 無限遠点を表す記号を文字 O から画像 に変更。Mon_Test と Mon_Test2 をわずかに改善。 (3) 追加の参考文献。 (4) §3の冒頭にイントロを追加。 (5) 表現の細部・表記の調整（約10カ所）。baby_steps テーブルの計算法について注釈 (*9) を追加。 (2) 「参考文献」: [Izu99] に注を付加。[Sut15] のフォーマットを他と統一。「日本語の文献」で「追加の参考文献」を紹介。各資料に西暦年数を付加。 (3) 表現の細部・表記の調整（約20カ所）。(*9) をサブセクション末尾に移動。P.length が1カ所 p.length となっていた。Ell_double の変数名 x1, y1 を x, y に変更。
