Hello.
ピピー！！👮👮無料版はここまでです！！！！！👮👮👮

---

# 521ビットのメルセンヌ素数を法とした剰余乗算
521ビットのメルセンヌ素数（$p = 2^{521}-1$）を法とする剰余乗算を効率的に行うお話です。
$GF(p):=\{0, 1, ..., p-1\}$とします。$x, y \in GF(p)$について$x$と$y$の積を$p$で割った余り、つまり$xy \,\,mod\,\, p$を計算したいのですが、剰余計算は乗算よりもコストが高いです。そこで、乗算と同じくらいのコストで剰余が求められるような方法を取り入れます。効率的な剰余のアルゴリズムと言うと、モンゴメリリダクション（以降MR）が有名だと思います。まずは一般的なMRについて紹介します。

法を$N$とし、$R>N$かつ$GCD(N, R)=1$となるような$R$を選択します。$N', R^{-1}$をそれぞれ$R^{-1} := 1/R \,\,mod\,\,N$, $N':= -1/N \,\,mod\,\, R$ とします。拡張ユークリッドのアルゴリズムを使うと$Ns + Rt =GCD(N, R)$を満たす$s,t$を得ることができ、$N'=-s\,\,mod\,\,R, R^{-1} = t \,\,mod\,\, N$です。これでMRの準備ができました。下に擬似コードを示します。

MRのアルゴリズム  
input: T(= xy)
output: TR^{-1} mod N

T <- (T + N(TN' mod R))/R
if T >= N then T <- T - N  
return T

MRを使うのは剰余計算を効率化するためでした。アルゴリズムを見ると$R$による除算と剰余計算が行われていてMRを使う意味がないように思われますが、$R$を2の冪にすることで剰余が論理積に、除算が論理右シフトになるので実は除算を行う必要がありません。出力が$xy\,\,mod\,\,N$ではなく$xyR^{-1}\,\,mod\,\,N$となっているため$xy\,\,mod\,\,N$を得るには$R$をかけてあげないといけないのですが、ここではMRのアルゴリズムのみの紹介に留め、これ以上は踏み込まないようにします。上の擬似コードで出力が正しく得られることの確認は皆さんに投げます。

さて、法が$p$（$=2^{521}-1$）の場合を考えてみましょう。$R$を$p+1=2^{521}$とします。$p$は素数なので$GCD(p,R)=1$を満たしています。$N', R^{-1}$を計算するとどちらも$1$です。これを反映させて擬似コードを再度示します。

MRのアルゴリズム(N=2^{521}-1の場合)  
input: T(= xy)
output: T mod p

T <- (T + p(T mod R))/R
if T >= p then T <- T - p  
return T

乗算が一回減り、出力は$xy \,\,mod \,\,p$になっています。
$(T\,\,mod\,\,R)$は$T$と$p$の論理積、$/R$は論理右シフトです。論理積を取ったあとの$p$倍は、任意の$r$について$rp = r(p+1)-r$が成り立ち$p+1=R$なので、$R$倍してから元の数を引いてあげればよいですね。$R$倍は$R$が$2$の冪なので論理左シフトです。もう一度、擬似コードを書きます。

MRのアルゴリズム(N=2^{521}-1の場合)  
input: T(= xy)
output: T mod p

T <- (T + R(T mod R) - (T mod R))/R
if T >= p then T <- T - p  
return T

剰余が足し算、引き算、論理シフトだけで計算できました。除算、乗算を使わずに剰余を求められるのは面白いですね。効率的に計算ができるのはコンピュータにとっても嬉しいです。

521ビットの楕円曲線（secp521r1）ではMRを利用して高速化を図ることができます。

## すこし応用
剰余乗算が効率的に行えるようになったので、次は冪剰余計算に挑戦したくなります。$p$が素数なので$p-2$乗してモジュラー逆数を求めてみましょう。

$p$はメルセンヌ素数ですので、$2$進数で書くと$(11..11)_{2}$です。ビットがすべて立っています。$p-2$は$(11..01)_{2}$になります。右向きバイナリ法が使えそうです。ビットパターンがわかっているので条件分岐をする必要がありません。擬似コードを書いてみます。

モジュラー逆数を求めるアルゴリズム
input: x
output: 1/x mod p

r <- x
for k in range(0, 518):
  r <- r*r mod p
  r <- r*x mod p

// 下位2ビット
r <- r*r mod p
r <- r*r mod p
r <- r*x mod p 
return r

乗算が519回、2乗が520回です。今後は乗算を$M$、2乗を$S$と書くことにします。計算コストは$519M+520S$となります。このコストを減らすためにアルゴリズムを改良していきましょう。シンプルな右向きバイナリ法だと$1$ビットずつ見るので、$k$ビットずつ見るようにすれば$519M$を$(1+\frac{518}{k})M$にできます。$k$ビットずつ見るアルゴリズムをsliding window methodとかk-aray methodと呼びます。本来のsliding window methodでは$x$の$0$乗から$2^{k}-1$乗までを事前計算しておく必要があるのですが、今回はビットがすべて$1$なので$x^{2^{k}-1}$だけを準備しておけばよさそうです。
519ビットを$k$ビットずつ見るということは$k$が519を割り切る必要があるように思えますが、$t:=519\,\,mod\,\,k$として初期値を$r:=x^{2^{t}-1}$とすることで柔軟に$k$を選択できます。(ただし$k>1$)
擬似コードを示します。$mod\,\,p$を省略してあります。

モジュラー逆数を求めるアルゴリズム（改良版）
input: x
output: 1/x mod p

r <- x^{2^{t}-1}
a <- x^{2^{k}-1}
for k in range(0, (519-t)/k):
  r <- r^{2^{k}}
  r <- r*a

r <- r*r 
r <- r*r 
r <- r*x 
return r

では計算コストを算出してみましょう。

事前計算: $(k-1)M+(k-1)S$
アルゴリズム本体: $(\frac{519-t}{k}+1)M+ 520S$

事前計算は工夫すればもっとコストを落とせるのですが、ひとまずこの式で算出することにします。アルゴリズム全体のコストが最小になるように$k$を探索すると、$k=17,18$付近で最小になることがわかります。

### 速度比較
（おそらく）最速なモジュラー逆数計算のアルゴリズムを紹介しました。
が、残念なことに各言語が提供しているpowm関数には勝てません。最速ではないことが判明しました。

C/C++を使う場合は、mpz_tならmpz_invertを、mp_limb_tならmpn_gcdextを使うのがよいと思います。さっき紹介したアルゴリズムを使ってもこの2つに速度で勝つのは相当難しいです。何度もmulModやsqrModを繰り返すので関数呼び出しのオーバーヘッドが大きいのだと思います（要因は他にもありそうですね）。

[実装と速度比較(ykm11 Gist)](https://gist.github.com/ykm11/c71a49ba92eb4e43cf4047d7d8bad6e7#file-pow521_v6-cpp)