\(p := 2^{521}-1 \): メルセンヌ素数
\(T := xy\), \(x, y \in F_{p} \)
\(T \, mod \, p\)を普通に計算してもいいけどせっかくなのでモンゴメリリダクションを使う.
今回はmodulusが特殊（全ビットが1）なので通常のモンゴメリリダクションよりも速くできる.
\(R := p+1 = 2^{521} \)
モンゴメリリダクションの式
\( (T + p(T \, mod \, R))/R = T \, mod \, p\)が成り立つ.
wikiとかには\(n^{'}\)や\(R^{-1}\)なるものが出てくるが今回はどちらも1なので書かない.
\(R\)は普通2のべきになるように取られる（modが論理積, 除算が右シフトになるので)が, modulusはプログラムによって変わるので戦略も変わる. 今回はかなり都合が良い.
\(T \, mod \, R\)のあとで\(p\)との乗算が発生しているが, \(p\)は全ビットが1であることと\(xp = x(p+1) - x\)が成り立つことを考えると左シフトと引き算だけで\(p\)との積を表すことができる.
足し算と引き算とシフト演算で剰余を求めることができるのは嬉しい. 爆速インターネットの誕生である.
速度比較はここ メルセンヌ素数2^{521}-1を法とする場合の乗算の実行時間比較 · GitHub
逆元の計算もバイナリ法でやってみた. \(p-2\)のビットパターンは既知なのでifを使わなくてもよい.
速度比較 メルセンヌ素数2^{521}-1の逆元計算の速度比較 · GitHub
mpn_invertには勝てなかった.
毎ビットでsqr, mulを呼び出すから, ここがオーバーヘッドになってるのかな.
上手いことwindowサイズ調整して乗算の回数減らしたら勝てそう.
secp521r1を実装するときはヤコビ座標で5bitsのNAFがいいかな（これより速いアルゴリズムもあるかもしれない）.
