C++でmod pの平方根を求める方法（Tonelli–Shanksアルゴリズム）

問題の概要

この問題では、整数 n と素数 p が与えられ、「n の mod p における平方根（モジュラ平方根）」を求めることを目標とします。

具体例を見てみましょう。

入力 : n = 4, p = 11
出力 : 9

これは、9² = 81 ≡ 4 (mod 11) が成立するためです。つまり 9 は 4 のモジュラ平方根となっています。

解法のアプローチ

本記事では Tonelli–Shanks（トネリ・シャンクス）アルゴリズムを使用します。

Tonelli–Shanksアルゴリズムとは、モジュラ演算において x² ≡ n (mod p) という形の合同式を満たす x を求めるための古典的かつ強力なアルゴリズムです。

アルゴリズムの手順

ステップ1 : n^((p−1)/2) mod p の値を計算します。結果が p − 1 と等しい場合、モジュラ平方根は存在しません（オイラーの基準による判定）。

ステップ2 : p − 1 を s × 2^e の形に分解します。ここで s は正の奇数、e は正の整数です。

ステップ3 : q^((p−1)/2) ≡ −1 (mod p) を満たすような非平方数 q を見つけます。

ステップ4 : ループ処理によって m を探索し、条件に応じて x などの変数を更新していきます。

まず、b^(2^m) ≡ 1 (mod p) となる m（0 ≤ m ≤ r−1）を探します。

m が 0 であれば x を答えとして返します。そうでない場合は、次のように各値を更新します。

x = x * g^(2^(r − m − 1))
b = b * g^(2^(r − m))
g = g^(2^(r − m − 1))
r = m

C++での実装例

上記の解法の動作を示すプログラムがこちらです。

#include <iostream>
#include <math.h>
using namespace std;
int powerMod(int base, int exponent, int modulus) {
   int result = 1;
   base = base % modulus;
   while (exponent > 0) {
      if (exponent % 2 == 1)
      result = (result * base)% modulus;
      exponent = exponent >> 1;
      base = (base * base) % modulus;
   }
   return result;
}
int gcd(int a, int b) {
   if (b == 0)
   return a;
   else
   return gcd(b, a % b);
}
int orderValues(int p, int b) {
   if (gcd(p, b) != 1) {
      return -1;
   }
   int k = 3;
   while (1) {
      if (powerMod(b, k, p) == 1)
      return k;
      k++;
   }
}
int findx2e(int x, int& e) {
   e = 0;
   while (x % 2 == 0) {
      x /= 2;
      e++;
   }
   return x;
}
int calcSquareRoot(int n, int p) {
   if (gcd(n, p) != 1) {
      return -1;
   }
   if (powerMod(n, (p - 1) / 2, p) == (p - 1)) {
      return -1;
   }
   int s, e;
   s = findx2e(p - 1, e);
   int q;
   for (q = 2; ; q++) {
      if (powerMod(q, (p - 1) / 2, p) == (p - 1))
      break;
   }
   int x = powerMod(n, (s + 1) / 2, p);
   int b = powerMod(n, s, p);
   int g = powerMod(q, s, p);
   int r = e;
   while (1) {
      int m;
      for (m = 0; m < r; m++) {
         if (orderValues(p, b) == -1)
         return -1;
         if (orderValues(p, b) == pow(2, m))
         break;
      }
      if (m == 0)
      return x;
      x = (x * powerMod(g, pow(2, r - m - 1), p)) % p;
      g = powerMod(g, pow(2, r - m), p);
      b = (b * g) % p;
      if (b == 1)
      return x;
      r = m;
   }
}
int main() {
   int n = 3;
   int p = 13;
   int sqrtVal = calcSquareRoot(n, p);
   if (sqrtVal == -1)
      cout<<"Modular square root is not exist";
   else
      cout<<"Modular square root of the number is "<<sqrtVal;
}

コードのポイント

• powerMod()： 繰り返し二乗法を用いて、大きな指数でも高速にべき剰余を計算する補助関数です。
• gcd()： 最大公約数を求めます。n と p が互いに素でない場合は解を返しません。
• orderValues()： mod p における b の位数（何乗すると 1 になるか）を求めます。
• findx2e()： 与えられた整数を奇数部分 s と 2 のべき乗部分 e に分解します。
• calcSquareRoot()： Tonelli–Shanksアルゴリズムの本体であり、これらの関数を組み合わせてモジュラ平方根を計算します。

出力

Modular square root of the number is 9

この実行結果から、9² ≡ 3 (mod 13) が成立していること、すなわち 9 が 3 のモジュラ平方根であることが確認できます。