牛牛喜欢看小姐姐[牛客挑战赛37F]

发表于 分类于

题目链接

牛客

题解

毒瘤数学题。

由裴蜀定理，不定方程 $ax+by=c$ 有整数解当且仅当 $\gcd(a,b)|c$。那么当一个小姐姐的步长为 $a_i$，能否走到 $x$ 的充要条件为 $\gcd(a_i,m)|x$。为了方便，下面的 $a_i$ 都变为 $\gcd(a_i,m)$。注意到 $a_i$ 都是 $m$ 的因子。而 $10^{18}$ 以内的因子个数是 $10^5$ 级别的。

考虑到只选 $s$ 个比较难统计。我们容斥一下，变成统计选至少 $s$ 个的答案。那么最后的答案就是至少 $s$ 个减去至少 $s+1$ 个。

那么现在就变成从 $n$ 个 $a_i$ 中选出 $s$ 个不同的数，记为 $b_i$。判断一个 $x$ 是否合法，当且仅当 $lcm(a_{b_1},a_{b_2}\cdots a_{b_s})|x$。我们需要统计出所有的满足条件的 $x$。对于这样一种方案，我们记集合： $A(lcm(a_{b_1},a_{b_2},\cdots,a_{b_s}))={lcm(a_{b_1},a_{b_2}\cdots a_{b_s})|x,x\in[1,k]}$。也就是所有满足条件的集合。

此时答案为：$\left | \bigcup_{b_1,b_2\cdots b_s}A(lcm(a_{b_1},a_{b_2}\cdots a_{b_s})) \right |$

集合的并集比较难处理。试试考虑用最普通的容斥原理，变成处理交集的情况？

一个点 $x$ 同时在一堆 $A(y_i)$ 中当且仅当 $lcm(y_i)|x$。

那么就有：$\left | \bigcap A(x_i) \right |=\left |A(lcm(x_i)) \right |$。因此交集的情况就好处理了。

只需要对于每个 $x$，算出容斥系数 $f_x$，答案就是 $\sum_{x|m}f_x|A_x|=\sum_{x|m}f_x\left \lfloor \frac{k}{x} \right \rfloor$。

这个 $x$ 实际上只有在 $x|m$ 的情况时才有意义，也就是 $10^5$ 级别。

现在的情况是，每次在 $\binom{n}{s}$ 个数中选出 $k$ 个数，这些数的lcm必须等于 $x$，贡献是 $(-1)^{k+1}$。

这个必须等于 $x$ 非常烦，考虑设 $g_x$ 为选出 $k$ 个数，这些数的lcm必须为 $x$ 的因数的答案。

发现这个 $f$ 的高维前缀和就是 $g$，那么算出 $g$ 之后高维差分一下就可以了。

显然，设 $y$ 为这 $\binom{n}{s}$ 个数中，是 $x$ 因数的个数。

那么有 $g_x=\sum_{i=1}^y\binom{y}{i}(-1)^{i+1}=(-1)(\sum_{i=0}^y\binom{y}{i}(-1)^i)+1=[y\not=0]$

考虑这个 $y$ 到底是什么，这 $\binom{n}{s}$ 个数是 $s$ 个数的lcm，这个 $lcm(a_{b_1},a_{b_2}\cdots a_{b_s})$ 要是 $x$ 的因数，也就是 $a_{b_i}$ 都要是 $x$ 的因数。显然当 $\sum_{i=1}^n[a_i|x]\geq s$ 时， $y\not =0$。

而这个 $\sum_{i=1}^n[a_i|x]$ 显然是一个高维前缀和的形式。

那么我们只需要用Pollard-Rho对 $m$ 进行质因数分解，暴力枚举出所有的因数，用哈希或者unordered_map 存起来。然后算高维前缀和以及高维差分就可以了。

程序

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
#include <bits/stdc++.h>
#include <unordered_map>
using namespace std;
#define FO(x) freopen(#x".in","r",stdin),freopen(#x".out","w",stdout)
#define fo(i,j,k) for(int i=(j),end_i=(k);i<=end_i;i++)
#define fd(i,j,k) for(int i=(j),end_i=(k);i>=end_i;i--)
#define DEBUG(x) cout<<#x<<"="<<x<<endl;
#define all(x) (x).begin(),(x).end()
#define cle(x) memset(x,0,sizeof(x))
#define ll long long
#define ull unsigned ll
#define db double
#define lb long db
#define pb push_back
#define mp make_pair
#define fi first
#define se second
inline ll read()
{
	ll x=0; char ch=getchar(); bool f=0;
	for(;ch<'0'||ch>'9';ch=getchar()) if(ch=='-') f=1;
	for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<3)+(x<<1)+(ch^48);
	return f?-x:x;
}

inline ll mul(ll a,ll b,const ll &p)
{
    ll ans=0;
    for(a%=p,b%=p;b;)
    {
        if(b&1) {ans+=a; if(ans>=p) ans-=p;}
        b>>=1; a+=a; if(a>=p) a-=p;
    }
    return ans;
}
inline ll Pow(ll a,ll b,const ll &p)
{
    ll ans=1;
    for(;b;b>>=1,a=mul(a,a,p)) if(b&1) ans=mul(ans,a,p);
    return ans;
}
ll gcd(ll x,ll y){return !y?x:gcd(y,x%y);}
ll ex_gcd(ll a,ll b,ll &x,ll &y)
{
    if(!b) {x=1; y=0; return a;}
    ll d=ex_gcd(b,a%b,y,x);
    y-=a/b*x; return d;
}
inline ll inv(ll a,ll m)
{
    ll x,y,d;
    d=ex_gcd(a,m,x,y);
    return (x+m)%m;
}
namespace MillerRabin{
	const int prime[10]={2,3,5,11,61,7,41,29,31,19};
	const int m=10;
	int t; ll r;
	inline bool witness(int a,ll n)
	{
		ll b=Pow(a,r,n);
		if(b==1) return 1;
		for(int i=0;i<t;i++,b=mul(b,b,n)) if(b==n-1) return 1;
		return 0;
	}
	inline bool isprime(ll n)
	{
		fo(i,0,m-1)
		{
			if(n==prime[i]) return 1;
			if(n%prime[i]==0) return 0;
		}
		r=n-1;
		for(t=0;!(r&1);r>>=1) t++;
		fo(i,0,m-1) if(!witness(prime[i],n)) return 0;
		fo(i,1,10) if(!witness(rand()%(n-1)+1,n)) return 0;
		return 1;
    }
}
using MillerRabin::isprime;
ll p[70]; int w[70]; int pn;
namespace PollardRho{
    ll rho(ll n)
	{
        ll c=rand()%(n-2)+2,x=rand()%n,y=x,d;
        for(int i=1,k=2;i++;)
        {
            x=(mul(x,x,n)+c)%n;
            d=gcd(abs(y-x),n);
            if(d!=1&&d!=n) return d;
            if(y==x) return n;
            if(i==k) y=x,k<<=1;
        }
    }
    void find(ll n)
	{
        if(n==1) return;
        if(isprime(n)) return (void)(p[++pn]=n);
        ll d=n; for(;d==n;) d=rho(n);
        find(d); find(n/d);
    }
}
using PollardRho::find;
const int N=3e5+5;
unordered_map<ll,int> ma;

ll d[N]; int dn;
ll n,m,k,s;
void dfs(int k,ll now)
{
	if(k>pn) {d[++dn]=now; return;}
	fo(i,0,w[k]) dfs(k+1,now),now*=p[k];
}
inline void get(ll m)
{
	find(m);
	sort(p+1,p+pn+1);
	int k=0;
	for(int i=1,j;i<=pn;i=j+1)
	{
		for(j=i;j<pn&&p[j+1]==p[i];j++);
		w[++k]=j-i+1; p[k]=p[i];
	}
	fo(i,k+1,pn) p[i]=0;
	pn=k;
	dfs(1,1);
	sort(d+1,d+dn+1);
	fo(i,1,dn) ma[d[i]]=i;
}
ll f[N],h[N];
inline ll solve(ll s)
{
	ll ans=0;
	int x=0;
	fo(i,1,dn) x=ma[d[i]],f[x]=(h[x]>=s);
	fo(i,1,pn)
		fd(j,dn,1)
			if(d[j]%p[i]==0)
				f[ma[d[j]]]-=f[ma[d[j]/p[i]]];
	fo(i,1,dn) ans+=f[ma[d[i]]]*(k/d[i]);
	return ans+(n>=s);
}

inline void calc_h()
{
	fo(i,1,n) h[ma[gcd(read(),m)]]++;
	fo(i,1,pn)
		fo(j,1,dn)
			if(d[j]%p[i]==0)
				h[ma[d[j]]]+=h[ma[d[j]/p[i]]];
}
int main()
{
	n=read(); m=read(); k=read(); s=read();
	get(m);
	calc_h();
	printf("%lld\n",solve(s)-solve(s+1));
	return 0;
}