https://ac.nowcoder.com/acm/contest/4381/F
先把(a[i])变成(gcd(a[i],m))，这样能到的就是a[i]的倍数。
应该第一眼想到容斥的，直接推系数太难了。
考虑容斥为至少s个，答案=至少s个-至少s+1个。
也就是选s个a[i]的lcm的倍数集合并起来的大小。
集合并可以用容斥搞成集合交。
一共有C(n,s)个lcm的倍数集合，选i个出来，取交，也就是所有lcm的lcm，容斥系数是((-1)^{i-1})。
不管怎样，lcm都是m的约数，所以枚举一个y|m作为最后的lcm。
设(f[y])为选若干个lcm的倍数集合，这些lcm的lcm=y的容斥系数和。
(g[y])表示这些lcm的lcm|y的容斥系数和。
(g)是(f)高维前缀和，(f)显然就是(g)的高维差分。
考虑求(g[y]) ，设(c=sum[a[i]|y])，(g[y]=sum_{i=1}^{C(c,s)}C(C(c,s),i)*(-1)^{i-1}=[C(c,s)>=1]=[c>=s])
(c)显然也可以高维前缀和求。
(Ans(至少s个)=sum_{d>=1~and~d|m}f[d]*(k/d)+[s<=n])
Code：
#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i < _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("
")
using namespace std;
const int N = 2e5 + 5;
ll n, m, k, s;
ll a[N];
ll gcd(ll x, ll y) {
return !y ? x : gcd(y, x % y);
}
ll mul(ll x, ll y, ll mo) {
x %= mo, y %= mo;
ll z = (long double) x * y / mo;
z = x * y - z * mo;
if(z < 0) z += mo; else if(z >= mo) z -= mo;
return z;
}
ll ksm(ll x, ll y, ll mo) {
ll s = 1;
for(; y; y /= 2, x = mul(x, x, mo))
if(y & 1) s = mul(s, x, mo);
return s;
}
int pd_p(ll n) {
if(n == 2) return 1;
if(n % 2 == 0) return 0;
ll s = n - 1, c = 0; while(s % 2 == 0) s /= 2, c ++;
fo(ii, 1, 40) {
ll x = ksm(rand() % (n - 1) + 1, s, n);
fo(i, 1, c) {
ll y = mul(x, x, n);
if(y == 1 && x != 1 && x != n - 1) return 0;
x = y;
}
if(x != 1) return 0;
}
return 1;
}
ll find(ll n) {
ll x = rand() % (n - 1) + 1, c = rand() % n, y = x;
ll i = 1, k = 2;
while(1) {
x = (mul(x, x, n) + c) % n;
ll d = gcd(n, abs(y - x));
if(d != 1 && d != n) return d;
if(x == y) return 1;
if((++ i) == k) y = x, k *= 2;
}
}
ll u[100], v[100], u0;
void fen(ll n) {
if(n == 1) return;
if(pd_p(n)) { u[++ u0] = n; return;}
ll d = find(n);
fen(d); fen(n / d);
}
void px() {
sort(u + 1, u + u0 + 1);
int U = u0; u0 = 0;
fo(i, 1, U) if(!u0 || u[i] != u[u0])
u[++ u0] = u[i], v[u0] = 1; else v[u0] ++;
}
ll d[200005]; int d0;
void dg(int x, ll y) {
if(x > u0) {
d[++ d0] = y;
return;
}
fo(i, 0, v[x]) {
dg(x + 1, y);
if(i < v[x]) y *= u[x];
}
}
const int M = 1960817;
ll h[M]; int h2[M];
int ha(ll n) {
int y = n % M;
while(h[y] != 0 && h[y] != n)
y = (y + 1) % M;
return y;
}
void add(ll n, int x) {
int y = ha(n);
h[y] = n; h2[y] = x;
}
int qu(ll n) {
return h2[ha(n)];
}
ll f[N], g[N];
#define ul unsigned long long
ul calc(ll s) {
ul ans = (s <= n ? 1 : 0);
fo(i, 1, d0) g[qu(d[i])] = (f[qu(d[i])] >= s);
fo(i, 1, u0) {
fd(j, d0, 1) if(d[j] % u[i] == 0)
g[qu(d[j])] -= g[qu(d[j] / u[i])];
}
fo(i, 1, d0) ans += (ul) g[qu(d[i])] * (k / d[i]);
return ans;
}
int main() {
srand(19260817);
scanf("%lld %lld %lld %lld", &n, &m, &k, &s);
fo(i, 1, n) scanf("%lld", &a[i]), a[i] = gcd(a[i], m);
if(m == 1) {
pp("%d
", s == n ? n : 0);
return 0;
}
u0 = 0; fen(m); px();
dg(1, 1);
sort(d + 1, d + d0 + 1);
fo(i, 1, d0) add(d[i], i);
fo(i, 1, n) f[qu(a[i])] ++;
fo(i, 1, u0) {
fo(j, 1, d0) if(d[j] % u[i] == 0)
f[qu(d[j])] += f[qu(d[j] / u[i])];
}
pp("%llu
", calc(s) - calc(s + 1));
}