题目
思路
考虑从a[1…n]枚举s个出来，求lcm，答案就是所有lcm倍数的并
一共有C(n,s)个lcm的倍数集合，选i个出来，取交，也就是所有lcm的lcm，容斥系数是(−1)i−1。
不管怎样，lcm都是m的约数
考虑枚举d|m，那么只需要算出lcm的lcm是它情况的容斥系数和就好了，记为f[d]
ans=f[d]*(k/d)
代码
#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("\n")
using namespace std;
const int N = 2e5 + 77;
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 yjy(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(time(NULL));
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\n",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\n",yjy(s) - yjy(s + 1));
}
