mod pでの平方根

QS での因数分解では \pmod{p} での平方根を利用します。つまりある整数 a (0\leqq a < p) について

を解きたいわけです。とはいえ 0 から p-1 まで全て 2 乗して確認する方法では明らかに効率が悪いです。この記事ではその計算を効率的に行う方法について記録しようと思います。(しばらく\pmod{p}の記載を省略します。)

平方剰余

そもそも剰余環での平方根は存在するかを確認しておく必要があります。というのも

なので 0\leqq x<p の p-1 個の数が 2 乗することで約半分の個数の数に対応する分けです。逆に言えば対応しない数もあるわけで、そのような数には平方根は存在しません。例えば \pmod{3} での 2 や \pmod{5} での 2,3 には平方根は存在しません。

ではその存在確認をどのようにするかというと、平方剰余というものを求めます。

という計算で、結果は 1 か -1 になります。(フェルマーの小定理 a^{p-1}\equiv 1 を移項して因数分解すると分かりやすいかもしれません)
この (a/p)=1 となる a を平方剰余といい、\pmod{p} で平方根が存在する数になります。一方 -1 となる a は非平方剰余といいます。

\bmod p での平方根

そして肝心の平方根の求め方ですが、以下のようになります。計算効率の観点から p \bmod 8 で場合分けします。

p \equiv 3, 7 \pmod8 の場合

1. x = a^{(p+1)/4} \bmod p
2. return x

p \equiv 5 \pmod8 の場合

1. x = a^{(p+3)/8} \bmod p
2. x^2 \not\equiv a \pmod p ならば
   1. x = x \cdot 2^{(p-1)/4} \bmod p
3. return x

p \equiv 1 \pmod8 の場合

1. d として非平方剰余を選ぶ
2. p-1=2^st として s, t を定める。ただし t は奇数
3. A \equiv a^t
4. D \equiv d^t
5. m = 0
6. For i in [0,s)
   1. if (AD^m)^{2^{s-1-i}} \equiv -1
      1. m = m + 2^i
7. x \equiv a^{(t+1)/2}D^{m/2}
8. return x

最後の p \bmod 8 = 1 でのアルゴリズムは Tonelli–Shanks のアルゴリズム (英語 Wikipedia) というアルゴリズムで、これ自体は他のケースにもそのまま適用することが可能です。また、最初の非平方剰余を選ぶ部分については[平方剰余]に書いた通り約半数の値が非平方剰余となるので 2 から 1 つずつ非平方剰余となるか調べて見つければすぐに見つかるはずです。(1 は必ず平方剰余なので飛ばす)
また、数式中に度々現れる剰余環でのべき乗については繰り返し二乗法で求められます。

int32 ModPow(int32 a0, int32 e, int32 p) {
  int64 r = 1;
  for (int64 a = a0; e; e >>= 1) {
    if (e & 1) r = r * a % p;
    a = a * a % p;
  }
  return r;
}

p \equiv 1 \pmod8 の場合その 2

書籍でもう 1 つ計算法が紹介されていたのでこちらも紹介します。まずアルゴリズムをそのまま書くと

1. t^2-a が非平方剰余となる t を選ぶ
2. x=(t+\sqrt{t^2-a})^{(p+1)/2} \pmod{p^2}
3. return x

となるわけですが、数式中に平方根が出てきています。これは\sqrt{t^2-a} の数値を出すということではなく \mathbf{F}_p[\sqrt{t^2-a}] で計算しましょうということです。\mathbf{F}_p[\sqrt{t^2-a}]？見慣れない式ですね。それもわかりにくいと思うのでとりあえず計算できる方法を書きましょう。
文字が多くなると大変なので \omega=\sqrt{t^2-a} と書きつつ整数 u_0, u_1, v_0, v_1 を使って積 (u_0+u_1\omega)(v_0+v_1\omega) を求めます。u,v などの計算は全て \bmod\ p してください。

何となく雰囲気に覚えはありませんか？複素数で i を導入したときに同様のことをしたと思いますし、\omega=i とすれば全く同じ計算になります。簡単に書くために (u_0,u_1)=u_0+u_1\omega とすると

となります。試していませんし証明できているわけでもないですが、恐らく (p+1)/2 乗することで \omega の係数が 0 になると思います。

template <typename T>
using IntInF = std::pair<T, T>;
// 動作未確認です
int32 PowModInF(IntInF<int32> a0, int32 e, int32 p) {
  IntInF<int64> r {1, 0};
  for (IntInF<int64> a {a0.first, a0.second}; e; e >>= 1) {
    if (e & 1) {
      r = IntInF<int64> {
        (r.first * a.first + r.second * a.second % p * omega) % p,
	(r.first * a.second + r.second * a.first) % p
      };
    }
    a = IntInF<int64> {
      (a.first * a.first + a.second * a.second % p * omega) % p,
      2 * a.first * a.second % p
    };
  }
  return r;
}

\bmod p^k での平方根

x^2\equiv a \pmod p ではなく x^2 \equiv a\pmod{p^k} の計算が必要になることがあるかもしれません。この計算も以下のように求めることができます。
まず x^2 \equiv a\pmod{p^k} を満たす x は x^2 \equiv a\pmod{p^{(k-1)}} も満たします。紛らわしいので \bmod\ p^k の解を x_k、\bmod\ p^{(k-1)} の解を x_{k-1} と書きましょう。すると適当な整数 n, m, l を使って以下のような関係式が成り立つことがわかります。

するとこの 3 式を混ぜることで

となります。この最初と最後を p^{(k-1)} で割ることで

さらに \bmod\ p の空間に放り込むことで

となり l が求まり、最初の 3 式の 3 つ目に代入することで x_k が求まります。

参考文献

• "Prime Numbers" by Richard Crandall (Springer, 2005) (邦訳版もあります)
• "Prime Numbers and Computer Methods for Factorization" by Hans Riesel (Birkhauser, 1994)