题解 P4935 【口袋里的纸飞机】
ComeIntoPower
·
·
题解
大家好我系出题人（还有个背锅侠mcfx）
题目背景是summer pockets!
对于测试点1,2，使用O(R^N\times N^2)的算法即可通过
对于测试点3,4,5,6，可以发现只需要枚举排序后的数列即可。然后使用一些组合数学计算答案。这样复杂度就变得可以接受。大约只有10^6种不同的数列
裸暴力的部分已经结束。现在我们来讲正解。考虑对每个数分开算贡献，即对于x\in[0,P-1]，算有x出现的网格数。
考虑补集转化一下，即算没有x出现的网格数。
x=0$可以特判，即如果没有$P$的倍数出现则不会出现$0
现在只考虑x\in[1,P-1]。没有x出现，相当于有若干对(s,t)(st\mod P=x)不能同时出现在数列里。由于P是质数，所以这样的有序对一定只有(P-1)对。
比如P=5,x=1，就有如下对：(1,1),(2,3),(3,2),(4,4)
这时，你发现除去s=t的，就相当于若干个互不干涉的限制：s和t不能同时出现。s=t的情况就是这个数本身不能出现。这意味着这些限制可以单独考虑方案，然后再全都做卷积
你又发现R\leq 10^9。这实在是太大了！不过你可以把余数相同的视为一体，假设这个余数出现了A次，就相当于如果选这个余数，就有A种方法。
可以发现所有余数的出现次数只有A和A+1这两种。所以刚才那些限制可以变成：
a0种(A,A)(即s,t都出现了A次),a1种(A,A+1),a2种(A+1,A),a3种(A+1,A+1)
所以，用指数生成函数表示，限制(A,A)的指数生成函数为e^{Ax}e^{Ax}-(e^{Ax}-1)(e^{Ax}-1)=e^{Ax}+e^{Ax}-1
答案的指数生成函数为(e^{Ax}+e^{Ax}-1)^{a0} (e^{Ax}+e^{(A+1)x}-1)^{a1+a2} (e^{(A+1)x}+e^{(A+1)x}-1)^{a3} e^{(R/P)x}（最后那个是考虑0的情况）
这已经可以快速计算，复杂度为O(n^2P+P^2)（P^2是处理a0,a1,a2,a3）
打表可以发现，本质不同的(a0,a1,a2,a3)只有O(\sqrt{P})种！这不好证明，但是在数据范围下是成立的。（mcfx：可以把这个序列当做是随机的，感性理解一下）
来自出题人的辣鸡标程（该标程可以大幅优化，但是懒得写了）
```cpp
// luogu-judger-enable-o2
#include<bits/stdc++.h>
#define maxP 5010
#define maxn 5010
using namespace std;
const int mod=1e9+7;
struct data{
int a,b,c;
data(){}
data(int a,int b,int c):a(a),b(b),c(c){}
bool operator<(const data& d)const{
if(a!=d.a)return a<d.a;
if(b!=d.b)return b<d.b;
return c<d.c;
}
};
map<data,int>st;
int sP,tg[maxP];
int ans,n,P,R,fac[maxn],inv[maxn],cnt[maxP][2];
int f[maxn],g[maxn],A,lim,tmp[maxn],dp0[250][maxn],dp1[250][maxn],dp2[250][maxn];
int sdp0[250][maxn],sdp1[250][maxn],sdp2[250][maxn];
void mul(int a[],int b[],int c[]){
for(int i=0;i<=n;++i)tmp[i]=0;
for(int i=0;i<=n;++i)
for(int j=0;i+j<=n;++j)
tmp[i+j]=(tmp[i+j]+1ll*a[i]*b[j])%mod;
for(int i=0;i<=n;++i)c[i]=tmp[i];
}
int qpow(int a,int b){
int ans=1;
for(;b;b>>=1,a=1ll*a*a%mod)
if(b&1)ans=1ll*ans*a%mod;
return ans;
}
int C(int x,int y){
if(y<0||x<y)return 0;
return 1ll*fac[x]*inv[y]%mod*inv[x-y]%mod;
}
void cp(int fr[],int to[]){
for(int i=0;i<=n;++i)to[i]=fr[i];
}
void get(int dp[250][maxn],int sdp[250][maxn],int a,int g[]){
cp(dp[a%sP],g);
mul(sdp[a/sP],g,g);
}
void init(int dp[250][maxn],int sdp[250][maxn]){
dp[0][0]=sdp[0][0]=1;
for(int i=2;i<=sP;++i)mul(dp[i-1],dp[1],dp[i]);
cp(dp[sP],sdp[1]);
for(int i=2;i<=sP;++i)mul(sdp[i-1],sdp[1],sdp[i]);
}
int main(){
scanf("%d%d%d",&n,&P,&R);
// fprintf(stderr,"[R%P=%d]",R%P);
fac[0]=inv[0]=1;
for(int i=1;i<=n;++i)fac[i]=1ll*fac[i-1]*i%mod;
inv[n]=qpow(fac[n],mod-2);
for(int i=n-1;i>=1;--i)inv[i]=1ll*inv[i+1]*(i+1)%mod;
A=R/P,lim=R%P;
for(int i=1;i<P;++i)
for(int j=i+1;j<P;++j){
// if(i*j%P!=2)continue;
if(i>lim&&j>lim)cnt[i*j%P][0]++; // 0: A
else if(i<=lim&&j<=lim)cnt[i*j%P][1]++; //1:A+1
}
for(int i=1;i<P;++i)tg[i*i%P]++;
for(int i=1;i<P;++i)
// if(i==2)
st[data(cnt[i][0],cnt[i][1],tg[i])]++;
sP=sqrt(P);
for(int i=0;i<=n;++i)dp0[1][i]=2ll*qpow(A,i)*inv[i]%mod;
dp0[1][0]--,init(dp0,sdp0);
for(int i=0;i<=n;++i)dp1[1][i]=1ll*(qpow(A,i)+qpow(A+1,i))*inv[i]%mod;
dp1[1][0]--,init(dp1,sdp1);
for(int i=0;i<=n;++i)dp2[1][i]=2ll*qpow(A+1,i)*inv[i]%mod;
dp2[1][0]--,init(dp2,sdp2);
for(map<data,int>::iterator it=st.begin();it!=st.end();++it){
int a0=it->first.a,a3=it->first.b,tg=it->first.c;
for(int i=0;i<=n;++i)f[i]=g[i]=0;
f[0]=1;
get(dp0,sdp0,a0,g),mul(f,g,f);
get(dp1,sdp1,(P-1-tg)/2-a0-a3,g),mul(f,g,f);
get(dp2,sdp2,a3,g),mul(f,g,f);
for(int i=0;i<=n;++i)g[i]=1ll*inv[i]*qpow(A,i)%mod;
mul(f,g,f);
ans=(ans-1ll*it->second*(f[n]))%mod;
}
ans=(1ll*fac[n]*ans+1ll*P*qpow(R,n)-qpow(R-A,n))%mod;
ans=(ans%mod+mod)%mod;
printf("%d",ans);
}
```