首先考虑对于一个 \(x\),什么样的数能够在 \(x\) 对应的集合中表示出来,不难发现一个数 \(y\) 属于 \(x\) 对应的集合,当且仅当其可以写成 \(x^{c_1b_1+c_2b_2+\cdots+c_mb_m}\) 的形式,而由于 \(p\) 是质数,根据费马小定理,指数上的值 \(\bmod(p-1)\) 后幂的值是不变的,因此根据斐蜀定理,记 \(d=\gcd(p-1,b_1,b_2,\cdots,b_m)\),那么对于某个 \(a_i\) 而言,其对应的集合就是 \(\{v|v=a_i^{kd},k\in\mathbb{Z}\}\)。
于是问题转化为 \(n\) 个形如 \(\{v|v=x^{kd},k\in\mathbb{Z}\}\) 的集合的并的大小。由于 \(p\) 是一个质数,我们考虑先找出 \(p\) 的一个原根 \(g\),然后求出每个数的离散对数 \(c_i\)——直接 BSGS 是 \(n\sqrt{p}\) 的无法通过,不过注意到我们 BSGS 的时间复杂度实际上是通过不等式 \(B+\dfrac{p}{B}\) 平衡出来的,故 \(n\) 次 BSGS 的复杂度则可以写成 \(B+\dfrac{np}{B}\) 的形式,使用均值不等式可知当 \(B=\sqrt{np}\) 时上式的值最小。
我们令 \(c_i=\gcd(c_i,p-1)\),问题又可以转化为,给定 \(n\) 个数 \(c_1,c_2,\cdots,c_n\),求 \(0\sim p-2\) 中有多少个数至少是一个 \(c_i\) 的倍数。预处理出 \(p-1\) 的全部因数,然后用类似于高维差分的方式容斥即可。
时间复杂度 \(p^{1/4}·\omega(p)·\log p+\sqrt{np}+d(p)^2\)。
const int MAXN = 1e5;
int n, m, mod, a[MAXN + 5], b[MAXN + 5], d, g, c[MAXN + 5];
int qpow(int x, int e) {
int ret = 1;
for (; e; e >>= 1, x = 1ll * x * x % mod)
if (e & 1) ret = 1ll * ret * x % mod;
return ret;
}
vector<int> pr;
bool check(int x) {
for (int p : pr) if (qpow(x, (mod - 1) / p) == 1) return 0;
return 1;
}
int hd[16777777], nxt[1111111], val[1111111], key[1111111], item_n = 0;
void ins(int x, int v) {
val[++item_n] = v; key[item_n] = x;
nxt[item_n] = hd[x & 16777215]; hd[x & 16777215] = item_n;
}
int query(int x) {
for (int e = hd[x & 16777215]; e; e = nxt[e]) if (key[e] == x) return val[e];
return -1;
}
int main() {
scanf("%d%d%d", &n, &m, &mod); d = mod - 1;
for (int i = 1; i <= n; i++) scanf("%d", &a[i]);
for (int i = 1; i <= m; i++) scanf("%d", &b[i]), d = __gcd(b[i], d);
for (int i = 1; i <= n; i++) a[i] = qpow(a[i], d);
int tmp = mod - 1;
for (int i = 2; i * i <= tmp; i++) if (tmp % i == 0) {
pr.pb(i);
while (tmp % i == 0) tmp /= i;
}
if (tmp > 1) pr.pb(tmp);
for (int i = 1; i <= min(500, mod - 1); i++) if (check(i)) {g = i; break;}
int pw = 1, lim = min(1000000, mod - 1);
for (int i = 0; i < lim; i++) ins(pw, i), pw = 1ll * pw * g % mod;
for (int i = 1; i <= n; i++) {
int ipw = qpow(pw, mod - 2), ppw = a[i];
for (int j = 0; j < 1000; j++) {
int p = query(ppw);
if (~p) {
c[i] = j * lim + p;
break;
}
ppw = 1ll * ppw * ipw % mod;
}
c[i] = __gcd(c[i], mod - 1);
}
vector<int> fac;
for (int i = 1; i * i <= mod - 1; i++) if ((mod - 1) % i == 0) {
fac.pb(i);
if ((mod - 1) / i != i) fac.pb((mod - 1) / i);
}
sort(fac.begin(), fac.end()); map<int, int> occ, cnt, has;
for (int i = 1; i <= n; i++) occ[c[i]] = 1;
for (pii p : occ) for (int j : fac) if (j % p.fi == 0) has[j] = 1;
int res = 0;
for (int j : fac) cnt[j] = (mod - 1) / j;
for (int i = (int)(fac.size()) - 1; ~i; i--) {
for (int j = i + 1; j < fac.size(); j++) if (fac[j] % fac[i] == 0) {
cnt[fac[i]] -= cnt[fac[j]];
}
}
for (pii p : has) res += cnt[p.fi];
printf("%d\n", res);
return 0;
}
/*
126 1 892371481
454064832 417570397 807650026 466892263 791344450 336924589 1944729 579965567 597532886 242884248 571201658 854252534 421251558 532865386 378268553 626243819 52655615 200538225 876249837 49382562 240030741 620587222 436393590 638009495 733238502 388289028 247259968 686579911 441669632 779962142 14713099 335151186 762271979 143469455 314598656 425836835 393234728 76558278 37681984 705123315 467409058 679552102 887234913 597098940 656080544 514409337 297942107 386863076 572681110 462982810 702355212 489867353 708204254 421873458 432787970 340529638 263956558 554893764 332152577 429576245 499254755 808916354 221523897 237603643 533279729 582300328 720434258 880524706 580775256 389508352 798343273 828515504 223160304 498268147 347717041 530130911 305412761 61678154 425081992 94690688 41647458 314689990 806098130 227769244 654240625 430235365 877551724 726687649 336637644 620085895 209835456 135659961 393357172 109025884 844121397 420061468 514991440 348027753 480616426 41993598 796546736 282041658 30961000 816644732 550127122 604405967 7475731 47513393 457741151 298634431 197223503 210065831 580653127 229891936 743203760 802277727 259142393 632894836 657867426 287851541 174610825 139623734 88981046 643076502 804691267 485717206
1
*/