luogu P5605 小 A 与两位神仙 - 原根 - 分治
题目传送门
Subtask 1
直接模拟。
Subtask 2
BSGS算法模板。
Subtask 3
考虑模 $m$ 的任意一个原根 $g$。
假设 $g^{ra} \equiv x \pmod {m}, g^{rb} \equiv y \pmod{m}$ 。
那么原题的方程等价于方程 $a \cdot ra \equiv rb \pmod {\varphi(m)}$。
它等价于 $x \cdot ra - y\cdot \varphi(m) = rb$。
设 $d = \varphi(m)$。
它存在满足条件的解当且仅当 $(ra, d) \mid rb$。
这个条件仍然不好处理，因为我们难以求出 $rb$。
考虑这样一个条件 $(ra, d) \mid (rb, d)$。
我们来证明，满足 $(ra, d) \mid rb$ 当且仅当 $(ra, d) \mid (rb, d)$。
充分性显然。
考虑必要性，因为 $((ra,d), (rb,d)) = (ra, rb, d) = ((ra, d) , rb) = (ra, d)$。
所以有 $(ra, d) \mid (rb, d)$。
现在的问题是如何求 $(ra, d)$。
引理 设 $g$ 为模 $m$ 意义下的一个原根 ，如果 $g^{ra} \equiv a \pmod {m}$，设 $\delta_m(a) = d$ ，那么有 $(ra, \varphi(m)) d = \varphi(m)$.
证明 因为有 $\left( g^{ra}\right)^d \equiv 1\pmod{m}$，所以 $\varphi(m) | ra \cdot d$。
那么有 $ra \cdot d \equiv 0 \pmod {\varphi(m)}$.
所以有 $d \equiv 0 \pmod{\frac{\varphi(m)}{(\varphi(m), ra)}}$.
因为 $d$ 是满足条件的最小正整数，所以有 $d = \frac{\varphi(m)}{(\varphi(m), ra)}$。
所以 $(\varphi(m), ra) d = \varphi(m)$。
因此我们把问题转化为求阶。
众所周知，阶是 $\varphi(m)$ 的约数，所以我们暴力枚举所有因子能通过这个subtask。
Subtask 4
因为 $p$ 随机，所以期望的因子个数只有几千，做法同 Subtask 3。
Subtask 5
骗暴力选手去卡常。
Subtask 6
我们考虑优化求阶的做法。
引理 设 $d = \delta_m(x)$，如果 $d_0$ 满足 $x^{d_0} \equiv 1\pmod {m}$，那么有 $d | d_0$。
证明 设 $d_0 = qd + r\ (0 < r < d)$。那么有 $x^{qd + r} \equiv (x^{d})^q x^r \equiv x^r \equiv 1\pmod{m}$。因为 $d$ 是最小的正整数满足 $x^{d} \equiv 1 \pmod {m}$，但 $r$ 是更小的满足条件的正整数，这和阶的定义矛盾。
初始令 $d =\varphi(m)$。根据欧拉定理，我们知道一定有 $x^{d} \equiv 1 \pmod{m}$。设 $d = k\delta_m(x)$
考虑枚举任意一个 $\varphi(m)$ 的质因子 $p$，如果满足 $p \mid d$，并且 $x^{d / p} \equiv 1 \pmod {m}$。 那么有 $p \mid k$，我们令 $d' = d / p$ 继续执行这个过程，直到不存在满足条件的 $p$。
由于乘法取模我们可以使用 __int128 ，所以时间复杂度 $O(T \log^2 m + n^{1/4})$。
如果同时求 $(ra, \varphi(m)), (rb, \varphi(m))$ 可能会被卡常。
Subtask 7
我们注意到，我们没有必要求 $(rb, \varphi(m))$。
设 $d = \varphi(m)$。
$(ra, d) | (rb, d) $ 等价于 $\frac{d}{(rb, d)}| \frac{d}{(ra, d)}$。
因此我们只用判断 $b^{e}$ 模 $m$ 的余数是否为1就行了。
其中 $e$ 是最小的正整数满足 $a^e \equiv 1 \pmod{m}$，即 $a$ 模 $m$ 的阶。
Subtask 7+
注意到这个 $O(\log^2 m)$ 求阶比较蠢，我们来优化一下。
设 $d = \prod_i p_i^{e_i}$。初始令 $d = \varphi(m)$。
考虑依次确定真实的 $e'_i$，对于每个 $p_i$，先让 $d' = \frac{d}{p_i^{e_i}}$，计算出 $v = x^{d}$，暴力枚举 $e'_i$，然后做 $v' = v^{p_i}$。
不难发现这个做法后面的暴力枚举的过程，快速幂均摊 $O(\log m)$。因此这个做法的复杂度是 $O(\omega(\varphi(m))\log m)$
注意到前面暴力枚举比较蠢，考虑套一个分治。定义一个函数 solve(l, r, y)，表示将区间左边的 $p^{e'}$ 和区间右边 $p^e$ 乘起来的结果记为 $t$，然后令 $y = x^t$，计算 $e'_l, \cdots, e'_r$。
考虑递归左边时把右边的 $p^e$ 乘起来做快速幂，递归右边的时候把左边的 $p^{e'}$ 乘起来做快速幂。
不难分析出这个做法的复杂度是 $O(\log \omega(\varphi(m)) \log m)$。
吐槽
- 不知道为什么很多人觉得 $(ra, d) \mid rb$ 当且仅当 $(ra, d) \mid (rb, d)$ 非常地显然，可能是我数学比较菜，觉得不是那么显然，需要说明一下。
- 成功区分了会求阶和不会求阶的选手
- 开始想卡掉求 2 次阶的做法，因为 F 很难，为了降低比赛难度，所以被要求这里不卡常，不过好像卡不卡区分度都是一样的。
- 事后发现，二进制分组和 long double 快速乘均可以轻松跑进两倍标程序时限内，感觉我对卡常一无所知，我似乎还忘了卡 long double 快速乘。（虽然我并不知道能个不能卡）。两倍常数还想卡？
- 洛谷怎么最大只支持 100 组数据
- 为什么当时众多验题人没人告诉我求阶还可以优化
Code
#include <bits/stdc++.h>
using namespace std;
typedef bool boolean;
#define ll long long
#define ull unsigned long long
template <typename T>
T add(T a, T b, T m) {
return ((a += b) >= m) ? (a - m) : (a);
}
template <typename T>
void pcopy(T* pst, const T* ped, T* pval) {
for ( ; pst != ped; *(pst++) = *(pval++));
}
ll mul(ll a, ll b, ll m) {
return ((__int128)a) * b % m;
/* ll rt = 0, pa = a;
for ( ; b; b >>= 1, pa = add(pa, pa, m))
if (b & 1)
rt = add(rt, pa, m);
return rt; */
}
ll qpow(ll a, ll p, ll m) {
ll rt = 1, pa = a;
for ( ; p; p >>= 1, pa = mul(pa, pa, m))
if (p & 1)
rt = mul(rt, pa, m);
return rt;
}
ll gcd(ll a, ll b) {
return (b) ? (gcd(b, a % b)) : (a);
}
ll randLL() {
static ull seed = 998244353, msk = (1ull << 61) - 1;
return (signed ll) ((seed = seed * seed + seed * 3 + 233) & msk);
}
boolean miller_rabin(ll n) {
static int T = 25;
static const int pri[10] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
if (!(n & 1))
return n == 2;
if (n < 1000) {
for (int p = 2; p * p <= n; p++)
if (!(n % p))
return false;
return true;
}
for (int i = 0; i < 10; i++) {
if (n == pri[i]) {
return true;
} else if (!(n % pri[i])) {
return false;
}
}
ll d = n - 1;
int s = 0;
while (!(d & 1))
s++, d >>= 1;
for (int t = 0; t < T; t++) {
ll b = randLL() % n;
if (!b)
continue;
ll tmp = qpow(b, d, n);
if (tmp == 1 || tmp == n - 1)
continue;
for (int i = 0; i < s; i++) {
tmp = mul(tmp, tmp, n);
if (tmp == n - 1)
goto nextTurn;
if (tmp == 1 || tmp == 0)
return false;
}
if (tmp != 1)
return false;
nextTurn:;
}
return true;
}
ll pollard_rho(ll x) {
// cerr << x << '\n';
ll a, b, c, g;
if (!(x & 1))
return 2;
while (true) {
b = a = randLL() % x;
c = randLL() % 127 % x;
do {
a = add(mul(a, a, x), c, x);
b = add(mul(b, b, x), c, x);
b = add(mul(b, b, x), c, x);
g = gcd(b - a, x);
(g < 0) ? (g = -g) : (0);
if (g == x)
break;
if (g > 1)
return g;
} while (a != b);
}
assert(false);
return 0;
}
void get_primary_factors(ll x, vector<ll>& rt) {
if (x == 1) {
return;
}
if (miller_rabin(x)) {
rt.push_back(x);
return;
}
ll a = pollard_rho(x);
get_primary_factors(a, rt);
get_primary_factors(x / a, rt);
}
vector< pair<ll, int> > get_primary_factor(vector<ll>& vec) {
vector< pair<ll, int> > rt;
if (vec.empty())
return rt;
sort(vec.begin(), vec.end());
vector<ll>::iterator it = vec.begin();
rt.push_back(make_pair(*it, 1));
for (it = it + 1 ; it != vec.end(); it++)
if (*it == rt.back().first)
rt.back().second++;
else
rt.push_back(make_pair(*it, 1));
return rt;
}
/// Template ends
typedef vector<pair<ll, int>> factor;
factor get_factor(ll m) {
vector<ll> tmp;
get_primary_factors(m, tmp);
return get_primary_factor(tmp);
}
int T;
ll m, phim;
factor fac;
class RankEvalutor {
public:
ll m;
vector<ll> vp;
vector<ll> pw;
vector<ll> pwd;
void init(ll m, factor& v) {
this->m = m;
int t = v.size();
vp.resize(t);
pw.resize(t);
pwd.resize(t);
for (int i = 0; i < t; i++) {
pw[i] = 1;
vp[i] = v[i].first;
for (int j = v[i].second; j; j--, pw[i] *= vp[i]);
}
}
void dividing(int l, int r, ll v) {
if (l == r) {
pwd[l] = 1;
while (v != 1) {
v = qpow(v, vp[l], m);
pwd[l] *= vp[l];
}
return;
}
int mid = (l + r) >> 1;
ll pd = 1;
for (int i = mid + 1; i <= r; i++) {
pd *= pw[i];
}
dividing(l, mid, qpow(v, pd, m));
pd = 1;
for (int i = l; i <= mid; i++) {
pd *= pwd[i];
}
dividing(mid + 1, r, qpow(v, pd, m));
}
ll eval(ll a) {
dividing(0, (signed) vp.size() - 1, a);
ll ret = 1;
for (auto x : pwd) {
ret *= x;
}
return ret;
}
} RankEvalutor;
int main() {
scanf("%lld%d", &m, &T);
fac = get_factor(m);
assert((signed) fac.size() == 1);
fac = get_factor(phim = m / fac[0].first * (fac[0].first - 1));
RankEvalutor.init(m, fac);
ll a, b, d;
while (T--) {
scanf("%lld%lld", &a, &b);
assert(gcd(a, m) == 1 && gcd(b, m) == 1);
assert(a < m && b < m && a >= 1 && b >= 1);
d = RankEvalutor.eval(a);
if (qpow(b, d, m) == 1) {
puts("Yes");
} else {
puts("No");
}
}
return 0;
}
浙公网安备 33010602011771号