bigint.js v0.5 beta 27:
JavaScriptでプチ楕円曲線暗号
点の数を数えて素数になるようにする（生成）

履歴

2004-08-28

48ビット、15桁

一週間を無駄にした初期設定ミスを直したら、 あっけなく、たった6万個台の b で成功、 JavaScriptだけを用いて、48ビットの楕円曲線（素数位数）が構成された。 コストは8時間半だった。ざっとチェックしたが間違いはないようだ。

法 281兆4749億7671万0731 は15桁だ。 実際には、上記の6万個のほか、別の範囲を1万個探しているが、 所要時間の8時間半は、その1万個も含んでいる。 ラッキーケースなのか、異様に速く見つかった。

y^2 = x^3 + 172358387239013*x + 200408220063094 (mod 281474976710731)
y = 11158390620848
Order = 281474976716851

2004-08-22

44ビット、14桁

44ビット（14桁）の生成に成功。 201素数を用い、見込み409554ステップ、期待値20万ステップのところ、12万7909ステップで成功した。11.6時間。

y^2 = x^3 + 2004082120250*x + 12345678127909 (mod 17592190240783) Order = 17592190238947

このサイズになると、最初の素数テーブルの準備に7分もかかった。 11.6時間はラッキーだったのか、それとも暫定経験式より実際は小さい期待値なのか。 最近のスケールでは、38ビット以降3例続けて期待値より小さめで成功している。

二進展開の共通部分をなるべく長くすることで高速化が期待できる。 0が多くなるようにする効果は、半分以上0にすることで、平均と比べて4:3程度の高速化、 半分以上1だけの場合と比べて、5:3程度の高速化と推定される。

48ビット、15桁

剰余演算で、割る数が14桁までだと、速いアルゴリズムがある。 48ビットは15桁になるため、剰余演算が遅くなる。 この原因によるコスト増加は、6～7割増しである。 ただし16桁以上にも適用できる一般的なアルゴリズムの場合よりは、半分程度で済む。

PigPGP 0.2.3 にも用いられている版の Bigint._mod_typeF に重大なバグがあった。修正した。

ECC.prototype._multiplyS2 の戻り値あたりで、this.add( comPoint, Result ) となっていたのは、 this.add2( comPoint, Result ) の間違い。 桁数が短いときはどちらでも同じだが、長いときはadd2にしないと結果が狂う。 この間違いにより計算精度が低下し、本来、無限遠点になって成功するはずのものが失敗していた可能性がある。 そうだとすれば全体のコストに悪影響を与えた。 36ビット以降について確認したが、現在の成功例については「誤差の影響で本当は失敗のものが成功と報告される」といったことはない。 つまり、たまたまそこが桁数の短い計算で無限遠点になった場合であろう。

2004-08-21

2004年8月21日 bigint.js v0.5 beta 27: JavaScriptでプチ楕円曲線暗号（更新）

点の位数を当てる現在のやり方では、もうあまり速くできない。 主にプラットフォームに依存する最適化で多少速くした。 12桁もでき、今日は13桁、40ビットまで行けた。

y^2 = x^3 + 865712345678*x + 333550052198 (mod 1099512345691)
mod は素数。解の数 1099512316009 も素数。両方1兆台（13桁）

これでもRSA256ビットくらいに相当するらしいので「プチ楕円曲線暗号」としては、 とりあえずいいかも。今のアルゴリズムでは、どっちにしろ生成は48ビットあたりが限界で60ビットはいけない。

覚え書き

引き続き、群の位数を素数にすることを続ける（現実の実装では素数でなくても、大きな素因数を持てばいい）。

現在のアルゴリズムでは、 b の捜索は分割して行える。 36ビットの法は、十分に現実的だ（数時間～運が良ければ10分で）。 数日がかりを覚悟するなら、48ビット程度までは、JavaScript で構成できると考えられる。

全体的に

点の数候補の素数のリストが長くなると、IEでは効率が落ちる。 平均的に、 候補リストがフルに近ければ、少しの b について調べるだけで成功するし、 候補リストを短くした場合、成功するまで、たくさんの b について調べなければならない。 しかし、1個の b、1個の素数あたりのコストで考えると、リストが大きめのとき 3ms なのに、小さめだと 2ms で約66%も時間が短くなる。 だから「リストを小さくして多くの b について調べる」方法は「リストを大きくして少しの b について調べる」方法より効率的だ。 新しい b になるごとに 30ms 程度のオーバーヘッドがかかるが、それを計算に入れても上記のとおりである。

法を長くしていけば、曲線の生成にかかる時間が伸びるのは当然であるが、 それは、ストレートには、点の数の候補の素数のリストが巨大化するためである。 とはいえ、ヘッセの不等式が示す範囲は法に比例するわけではなく、 さいわい0.5乗（平方根）がくっついているので、法を長くしても急激には困難にならない。 さらに、数が大きくなれば素数自体の分布密度も少しずつ薄まるため、点の数の候補の増大はゆるやかである。 しかしながら、プラットフォームによっては、 リストが大きくなると、大きなリストを保持していることそれ自体が速度に悪影響を及ぼす。 1回あたりの発見確率を犠牲にしても小さなリストにして、回数を増やすほうが得策だ。

もうひとつは、法が伸びれば、当然、ビット単位でのかけ算や足し算のコストが増加する。 影響が大きい繰り返し二乗法の部分については、二進展開したときに1が少なくなるように、点の数の候補を絞っておくことで、 多少、改善が期待される。

成功するまでに試さなければいけない回数の期待値

経験的に、ごくおおざっぱではあるが、下記の hope 個の b （平方非剰余でスキップされる分もカウント）について調べれば、だいたい1回は成功する。
var hope = Math.log( modulus ) * 7 * ( 1.4 / hs ) * ( len / idx );
ここで、hs はヘッセの不等式の係数（本来は2）をどこまで絞るかである。 絞れば絞るほど、調べるべき b は増大する。 ヘッセの不等式が示す範囲の両端にはあまり分布していないこともあり、2に近い数にすると効率が落ちる。 len はそうして決定された範囲に存在する素数の個数である。 idx は len 個ある素数のうち、実際に何個を使用するかである。 実装では、高速化を希望して、二進展開したとき1が多い素数は捨てる、という作業を行うため、idx は len よりかなり小さい。

期待値としては、上記の半分の値が一応の目安となる。

11桁、36ビット

法を36ビットの素数とする素数位数の楕円曲線の構成に、2回成功している。36ビットの曲線生成は、JavaScriptでも十分にコントロール可能と考えられる。

最初の試行では、 hs を 0.4 とし、範囲内の8369個の素数のうち、 2449個を使用した。 上記の経験的評価式によると、 3588個の b についてチェックする必要がある。 その半数を平方剰余と考えて良い。 実測によると1個あたり7.4秒であり、 約3.7時間（期待値1.9時間）を見込まなければならない。実際には4.5時間かかり、 約4300個の b を試している。だいたい予測時間のとおりであった。

効率についての試行錯誤のあと、 hs を0.14まで絞り、範囲内2944個の素数のうち、わずか121個を使用するようにした。 これだと、42499個もの b について調べなければならないが、実測では、平方剰余の b 1個につき、 わずか230.5msであるから、上記と同様の考察によれば、1.4時間を見込めば良いことになる。 半分以下の時間でできることが期待される。 実際には4180個の b でけりがつき、わずか506秒であった。 ラッキーケースで、見込みの10分の1の時間で良い b に遭遇した。

12桁、38ビット

係数を0.1まで絞り、3954個の素数中、 238個だけを使う。 上記経験式によると、成功までに42884個の b について試す必要がある。 実測によると、平方剰余の b 一つが551msであり、 3.3時間が見込まれる（期待値1.7時間）。実際には、20516個のbを試したところで見つかった。約1.6時間であった。 

y^2 = x^3 + 123456789012*x + 90981020515 (mod 274977907007)
y = 58752894331
Order = 274977858061

13桁、40ビット

係数を0.05に絞り、 107 / 3783個の素数を使用。 19万2129個の b を試すうちには成功が見込まれる。 約6時間の見込みだが（期待値3時間）、実際には5万2198個、1.7時間で成功した。 

**** Prime Order ***
y^2 = x^3 + 865712345678*x + 333550052198 (mod 1099512345691)
y = 214322749639
Order = 1099512316009

2004-08-20

法9桁に成功。

y^2 = x^3 + 123456789*x + 89012095 (mod 234567899)

y^2 = x^3 + 123456789*x + 321099074 (mod 505059089)

下記は、法10桁。1ステップ15秒だった。

y^2 = x^3 + 876543210*x + 505051080 (mod 1234567891)

下記は、法10桁、32ビット。ステップ25.5秒。所要17分。

y^2 = x^3 + 1720304083*x + 2718281826 (mod 4321007869)

下記は運が良かった場合で所要5分。

y^2 = x^3 + 1720304055*x + 2718281886 (mod 4321007869)

33ビット。 ……点の数の候補の素数を二進展開したときの1の数は、 繰り返し二乗法でのかけ算の部分の計算量に相当する。 もしハッセの不等式の範囲を全数検索せず運任せで範囲を絞るなら、 二進展開で1の数が多いものはどんどん捨てたほうがいい。 発見までの平均ステップ数が伸びる代わり、1ステップは非常に速くなる。 発見までに必要な平均総ステップ数はどちらにしても変わらないが、かけ算の便利さと無関係に、点の数は均質に分布しているから、 速く探せるところばかり探せば結局コストダウンになる。 まじめにやった場合と比べて、0.6倍程度の時間でできるようだ。

y^2 = x^3 + 6400220064*x + 3031991588 (mod 9999999943)

11桁、34ビット。1ステップ10秒で約150ステップで見つかった（30分）。

y^2 = x^3 + 4321043210*x + 15697988204 (mod 17180868043)

36-bit。「手抜き法」を使ったため約4300のbを検査することになった。 発見までに実質4.5時間かかっている。1ステップ7400ms。 4000上から始めていたら300で見つかったのだが…

y^2 = x^3 + 27182818284*x + 31415930324 (mod 68732103211)

2004-08-19

このデモは、素数位数の有限体上に楕円曲線を作って、その曲線の位数も素数になるように曲線の定数項を決める。 現在、法が8桁まで実用になる。

「JavaScript上での楕円曲線暗号について」（2003年12月8日）の改善。 一番影響が大きいのは次のトリック。 点にかけるスカラーの最小と最大が二進展開して
11101101001110110000001
11101101011111000010101
とする。上位約40%の
11101101000000000000000
は、どのスカラー倍でも共通であるから、これを「くくり出して」使い回すことで約40%高速化する。

高速化 他の整理などと合わせて、従来の約2倍で生成できる。 法が7桁のとき、1つのb（平方剰余とする）を試す（そして曲線の位数が合成数と確認する）のに、約650ms。 50ステップで素数位数の曲線が見つかるとすると、30秒ほど。

8桁の法 上記アルゴリズムで、支障なく法8桁に進める。 JavaScript自身の演算で掛けてもオーバーフローしない26.5ビットまで容易に到達できる。 JavaScript自身の限界ぎりぎりの
y^2 = x^3 + 87654321*x + 54545454 (mod 94906247)
を構成する場合、上記の1ステップが約2.5秒。

9桁の法 少なくとも、積をとる部分は自力でオーバーフローを回避しなければならない。 原理的には演算を全部 Bigint でやれば問題ないが、それだとかなり遅い。 積以外の和などは（一応15桁まで）オーバーフローしないので、15桁程度までは、なるべく組み込みの足し算などを使うことで、 速度の低下を抑えられるはずだ。現在、うまく実装できていない。

メモ

Bigint.getPrimeTable( 7777000, 7777999 ) で素数表を作り、4で割ると3余る素数をどれか一つ選ぶ。 ここでは 7777727 とした。これが法である。次に x の係数を勝手に選び固定する。ここでは 2345678 とした。 この条件で定数項を変化させて、点の数が適切になるようにする。具体的には
y^2 = x^3 + 2345678x + b ( mod 7777727 )
の上に乗る点の数（位数）を数える。

ここでは位数を素数にしたい。 素数であるとすると、Hasse(-Weil)の定理によって、その素数は
(7777727 + 1) ± 2(7777727)^(1/2) の範囲にある。具体的に、次のようにして、715要素からなる素数表を手に入れる。 

    var modulus = 7777727;
    var min = Math.ceil( ( modulus + 1 ) - 2 * Math.sqrt( modulus ) );
    var max = Math.floor( ( modulus + 1 ) + 2 * Math.sqrt( modulus ) );
    var myPrimes = Bigint.getPrimeTable( min, max );

b は法より小さければ何でもいいが、ここでは b = 3000000 から始めて、良い結果が出るまで 1 ずつ増やして試行しよう。 ヤコビ記号で平方剰余か調べて、そうならば平方根を求める。そうでなければ次へ進む。 平方根が求まった場合、つまり曲線上の非自明な点がひとつ分かった場合、その点に位数候補の素数を次々とかける。 ここでほぼ同じ大きさのスカラーを繰り返し同じ点にかけるのだから、 繰り返し二乗法の計算の相当の部分が同じ反復になるのであって、 途中結果をキャッシュすることである程度の高速化が期待される。 それは後で考えることにして、ここではベタでやっておく。 すなわち、まずは、2003年12月 8日のメモにある「g_multiplyS_Work テーブル」を使う。 そのため、myPrimesテーブルの最初と最後の要素だけスワップして、最初の要素がいちばん大きくなるようにしておく。 

    var minPrime = myPrimes[ 0 ];
    var maxPrime = myPrimes[ len - 1 ];
    myPrimes[ 0 ] = maxPrime;
    myPrimes[ len - 1 ] = minPrime;

50秒で b = 3000054 素数位数 7778291 が見つかった。 さらに続けると、48秒で b = 3000118 : 素数位数 = 7780637 が、 68秒で 3000203 : 7773761 が、 92秒で 3000310 : 7778149 が、 45秒で 3000362 : 7773187 が、 142秒で 3000557 : 7780561 が、 94秒で 3000673 : 7774399 が、それぞれれ見つかった。

だいたいこのスケールでは、b を100個試すうちに1個くらいは素数位数になる。 ならない場合、1個の b を試すのに1550ms程度かかった。 高速化のために、まず二進文字列を無意味にひっくり返していた場所があるので、 ひっくり返さないではじめから逆さまに読むようにしたところ、 5%強、速くなり、1450ms にあがった。

次に、b を変化させながら、その累乗をいろいろ計算するするが、累乗の指数は myPrimes リストを繰り返し使う。 つまり、そのたびに同じ二進展開をするのは無駄なので、myPrimesBinString テーブルを一回だけ作って、あとはそれを使い回すほうが、 計算ステップは減る。しかし、myPrimesBinString テーブルはやや大きくなるので、あまり速度向上には結びつかず、わずかに逆効果だった。

ハッセの不等式の係数を 2 から1.5に絞る。素数位数の検出精度はわずかに下がるが、非素数位数が素数と誤認されることはない。 探す範囲が75%になったので、25% 高速化して、1段階が1100msになった。 係数を 1 にすれば当然50%高速化して、一段階700msになる。 それだと 7774957 から 7780511 までしか探さないから、最初にある素数位数の7つのうち、2つしかヒットしないことになって、 見逃す率が高すぎる。係数が1.5のときは、 7773569 から 7781909 まで探し、 7例のうち6つはヒットする。 つまり係数を1にすると速度的にはさらに25%ほど有利だが、その代償として発見までのステップが倍以上になる可能性が高いから、 速度面で欲張りすぎて損をする形だ。係数は1.5～1.8程度が安全であろう。 ここでは1.5とする。

すると、点にかけるスカラーの最小と最大は、二進展開して
11101101001110110000001
11101101011111000010101
となる。 上位約40%の
11101101000000000000000
は、どのスカラー倍でも共通であるから、これを「くくり出して」使い回すことで、節約が考えられる。

すなわち、次のようにして、共通因子をとりあえずくくりだし、 Point( 0, y ) に大きなスカラー p をかける代わりに、 最初に共通因子をかけて共通の点 ComPoint を得ておき、 p から共通因子を引いた小さい因子をPointにかけ、結果の点と ComPoint をかければ、 p をかけたことになる。 

    var comBinString = "";
    for( var pos = 0; pos < strlen; pos++ )
    {
        var c = maxBinString.charAt( pos );
        if( minBinString.charAt( pos ) !== c ) break;
        else comBinString += c;
    }

    for( i = pos; i < strlen; i++ ) comBinString += "0";

例えば、Point に二進数で
11101101001110110000001
11101101011111000010101
などをかける代わりに、最初に Point に共通部分の
11101101000000000000000
をかけて ComPoint を得ておき、ComPoint に非共通部分の 001110110000001
011111000010101
などをかけることで、目的を達する。

このトリックによって、 必要なかけ算の量が半分近くになり、40%近く高速化し、 1ステップが1100→700msになった。 このトリックが簡単に働くためには、最小の素数と最大の素数がそれぞれ二進展開して同じ長さになる必要がある。 万一そうなっていないときは、法を変更する。 また、非共通部分を文字列の配列で保持するのは、前述のように実装上必ずしも効率が良くない面もあるが、 どっちにしても二進展開した文字列が後で必要になること、 そして二進展開の文字列比較の代わりに2で割り算を繰り返して数値的に「共通因子」の検出をする方法もそれなりにコストがかかること、 などから、これはこのままとする。（このスケールでは、 文字列配列の作成時間は30ミリ秒であり、とりあえず良いとする。）

非共通部分を格納した myBinString リストでは文字列が0から始まる場合もある。 不要な0だが、除去するのは面倒であり、あってもほとんど悪影響がないので、そのままとする。 非共通部分はすべて同じ長さであるから、 myBinString がどういう順序の表でも「g_multiplyS_Work テーブル」は正しく作成される。 したがって、myPrimesリストの最初の要素が最大になるようなスワップも、もはや必要なくなる。 

    for( i = 0; i < len; i++ )
    {
        myBinString[ i ] = ( myPrimes[ i ].toString( 2 ) ).substr( pos );
    }

以上の変更によって b = 3000000 からインクリメントして素数位数を発生させる 3000054 を検出するまでのコストは、 22.5秒となった。最初は50秒かかっていたので、半分の時間で動くようになった。

Demos1
Demos2
Demos 3 (Slow)
Demos 4 (Slow)
Demos 5
Demos 6
Output
Debug
faireal.net

[index]