cp-library

C++ Library for Competitive Programming

View the Project on GitHub emthrm/cp-library

:heavy_check_mark: 平方剰余 (quadratic residue)
(include/emthrm/math/mod_sqrt.hpp)

平方剰余 (quadratic residue)

$x^2 \equiv a \pmod{p}$ を満たす $x$ が存在すれば、$a \in \mathbb{Z}$ は法 $p$ の下で平方剰余であり、そうでなければ平方非剰余である。

ルジャンドル記号 (Legendre symbol)

整数 $a$、奇素数 $p$ に対して

\[\left(\dfrac{a}{p} \right) \mathrel{:=} \begin{cases} 1 & (a \not\equiv 0 \pmod{p} \wedge a \text{ は法 } p \text{ の下で平方剰余}), \\ -1 & (a \text{ は法 } p \text{ の下で平方非剰余}), \\ 0 & (a \equiv 0 \pmod{p}) \end{cases}\]

と定義する。

オイラーの基準 (Euler’s criterion)

整数 $a \neq 0$、奇素数 $p$ に対して $a \perp p$ ならば

\[\left(\dfrac{a}{p} \right) \equiv a^{\frac{p - 1}{2}} \pmod{p}\]

が成り立つ。

Tonelli–Shanks algorithm

整数 $a$、奇素数 $p$ に対して $x^2 \equiv a \pmod{p}$ が成り立つ $x$ ($0 \leq x < p$) を求めるアルゴリズムである。

以下では剰余演算のときに $p$ を法とする。

  1. 一般性を失わず $0 \leq a < p$ とできる。

    $a = 0$ ならば $x = 0$ として終了する。

    $a$ が平方非剰余ならば条件を満たす $x$ は存在しないとして終了する。

  2. $p = q2^s + 1$ を満たす奇数 $q$、正整数 $s$ を求める。

    $p$ は奇素数より $q, s$ は一意に決まる。

  3. 平方非剰余 $z$ を一つ求め、$m = s,\ c \equiv z^q,\ t \equiv a^q,\ r \equiv a^{\frac{q + 1}{2}}$ とする。このとき

    • $c^{2^{m - 1}} \equiv -1 \quad (\because z \text{ は平方非剰余} \wedge c^{2^{m - 1}} \equiv z^{q \cdot 2^{m - 1}} = z^{\frac{p - 1}{2}})$、

    • $t^{2^{m - 1}} \equiv 1 \quad (\because a \text{ は平方剰余} \wedge t^{2^{m - 1}} \equiv a^{q \cdot 2^{m - 1}} = a^{\frac{p - 1}{2}})$、

    • $r^2 \equiv at \quad (\because r^2 \equiv a^{q + 1})$

    が成り立つ。この3式を成り立たせつつ $t \equiv 1$ となる $m, c, r$ を求めればよい。

  4. $t \equiv 1$ ならば $x = r$ として終了する。

  5. $t^{2^i} \equiv 1$ を満たす最小の正整数 $i$ を求める。

    上の2式目より $i = m - 1$ はこれを満たすので $1 \leq i \leq m - 1$ を満たす。

    $t \not\equiv 1$ より $t^{2^{i - 1}} \equiv -1$ が成り立つ。

    • $m \leftarrow i$、

    • $c \leftarrow c^{2^{m - i}} \quad (\because (c^{2^{m - i}})^{2^{i - 1}} = c^{2^{m - 1}} \equiv -1)$、

    • $t \leftarrow tc^{2^{m - i}} \quad (\because (tc^{2^{m - i}})^{2^{i - 1}} = t^{2^{i - 1}}c^{2^{m - 1}} \equiv -1 \cdot -1 \equiv 1)$、

    • $r \leftarrow rc^{2^{m - 1 - i}} \quad (\because (rc^{2^{m - 1 - i}})^2 \equiv r^2 c^{2^{m - i}} \equiv atc^{2^{m - i}})$

    と更新しても3式は成り立つ。この操作で $m$ は減少している。$m = 1$ のとき2式目より $t \equiv 1$ が成り立つため、有限ステップで解は求まる。

  6. 4に戻る。

ヤコビ記号 (Jacobi symbol)

整数 $a$、正の奇数 $p$ に対して $p$ の素因数分解を $p = \prod_i p_i^{e_i}$ とすると

\[\left(\dfrac{a}{p} \right) \mathrel{:=} \prod_i \left(\dfrac{a}{p_i} \right)^{e_i}\]

と定義される。

時間計算量

  時間計算量
平方剰余  
ヤコビ記号 $O(\log{\max \lbrace A, P \rbrace})$ ?

仕様

平方剰余

名前 戻り値 要件
long long mod_sqrt(long long a, const int p); $x^2 \equiv a \pmod{p}$ を満たす $x$。ただし存在しないときは $-1$ を返す。 $p \in \mathbb{P}$

ヤコビ記号

名前 戻り値
int jacobi_symbol(long long a, long long p); $\left(\dfrac{a}{p} \right)$

参考文献

平方剰余

ルジャンドル記号

ヤコビ記号

TODO

Submissons

Depends on

Verified with

Code

#ifndef EMTHRM_MATH_MOD_SQRT_HPP_
#define EMTHRM_MATH_MOD_SQRT_HPP_

#include <random>

#include "emthrm/math/mod_pow.hpp"

namespace emthrm {

long long mod_sqrt(long long a, const int p) {
  if ((a %= p) < 0) a += p;
  if (a == 0) [[unlikely]] return 0;
  if (p == 2) [[unlikely]] return 1;
  if (mod_pow(a, (p - 1) >> 1, p) == p - 1) return -1;
  if (p % 4 == 3) return mod_pow(a, (p + 1) >> 2, p);
  int s = 1, q = (p - 1) >> 1;
  for (; !(q & 1); q >>= 1) {
    ++s;
  }
  static std::mt19937_64 engine(std::random_device {} ());
  std::uniform_int_distribution<> dist(2, p - 1);
  long long z;
  do {
    z = dist(engine);
  } while (mod_pow(z, (p - 1) >> 1, p) == 1);
  int m = s;
  long long c = mod_pow(z, q, p), r = mod_pow(a, (q - 1) >> 1, p);
  long long t = a * r % p * r % p;
  r = (r * a) % p;
  while (t != 1) {
    long long t2 = t * t % p;
    for (int i = 1; i < m; ++i) {
      if (t2 == 1) {
        const long long b = mod_pow(c, 1 << (m - i - 1), p);
        m = i;
        r = (r * b) % p;
        c = b * b % p;
        t = (t * c) % p;
        break;
      }
      t2 = (t2 * t2) % p;
    }
  }
  return r;
}

}  // namespace emthrm

#endif  // EMTHRM_MATH_MOD_SQRT_HPP_
#line 1 "include/emthrm/math/mod_sqrt.hpp"



#include <random>

#line 1 "include/emthrm/math/mod_pow.hpp"



namespace emthrm {

long long mod_pow(long long x, long long n, const int m) {
  if ((x %= m) < 0) x += m;
  long long res = 1;
  for (; n > 0; n >>= 1) {
    if (n & 1) res = (res * x) % m;
    x = (x * x) % m;
  }
  return res;
}

}  // namespace emthrm


#line 7 "include/emthrm/math/mod_sqrt.hpp"

namespace emthrm {

long long mod_sqrt(long long a, const int p) {
  if ((a %= p) < 0) a += p;
  if (a == 0) [[unlikely]] return 0;
  if (p == 2) [[unlikely]] return 1;
  if (mod_pow(a, (p - 1) >> 1, p) == p - 1) return -1;
  if (p % 4 == 3) return mod_pow(a, (p + 1) >> 2, p);
  int s = 1, q = (p - 1) >> 1;
  for (; !(q & 1); q >>= 1) {
    ++s;
  }
  static std::mt19937_64 engine(std::random_device {} ());
  std::uniform_int_distribution<> dist(2, p - 1);
  long long z;
  do {
    z = dist(engine);
  } while (mod_pow(z, (p - 1) >> 1, p) == 1);
  int m = s;
  long long c = mod_pow(z, q, p), r = mod_pow(a, (q - 1) >> 1, p);
  long long t = a * r % p * r % p;
  r = (r * a) % p;
  while (t != 1) {
    long long t2 = t * t % p;
    for (int i = 1; i < m; ++i) {
      if (t2 == 1) {
        const long long b = mod_pow(c, 1 << (m - i - 1), p);
        m = i;
        r = (r * b) % p;
        c = b * b % p;
        t = (t * c) % p;
        break;
      }
      t2 = (t2 * t2) % p;
    }
  }
  return r;
}

}  // namespace emthrm
Back to top page