C++ Library for Competitive Programming
#include "emthrm/math/mod_sqrt.hpp"
$x^2 \equiv a \pmod{p}$ を満たす $x$ が存在すれば、$a \in \mathbb{Z}$ は法 $p$ の下で平方剰余であり、そうでなければ平方非剰余である。
整数 $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}\]と定義する。
整数 $a \neq 0$、奇素数 $p$ に対して $a \perp p$ ならば
\[\left(\dfrac{a}{p} \right) \equiv a^{\frac{p - 1}{2}} \pmod{p}\]が成り立つ。
整数 $a$、奇素数 $p$ に対して $x^2 \equiv a \pmod{p}$ が成り立つ $x$ ($0 \leq x < p$) を求めるアルゴリズムである。
以下では剰余演算のときに $p$ を法とする。
一般性を失わず $0 \leq a < p$ とできる。
$a = 0$ ならば $x = 0$ として終了する。
$a$ が平方非剰余ならば条件を満たす $x$ は存在しないとして終了する。
$p = q2^s + 1$ を満たす奇数 $q$、正整数 $s$ を求める。
$p$ は奇素数より $q, s$ は一意に決まる。
平方非剰余 $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$ を求めればよい。
$t \equiv 1$ ならば $x = r$ として終了する。
$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$ が成り立つため、有限ステップで解は求まる。
4に戻る。
整数 $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)$ |
平方剰余
ルジャンドル記号
ヤコビ記号
#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