牛客挑战赛37 F 牛牛喜欢看小姐姐（容斥+高位前缀和）

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));
}