「USACO 2020 US Open Platinum」Exercise
做法与模数是否是质数无关
问题可能比较复杂，需要多步分析
1.对于一个已知的排列
显然这样的置换会构成若干个环，设每个环长度为(a_i,iin [1,m])，显然答案就是(lcm(a_i))
2.对于已知的(a_i)序列(注意这里说的是有序的)，计算其方案数
考虑已经排列的个数为(i)，加入一个环大小为(j)
为了避免重复，应当固定这个环的初始位置为1号点，其余位置按照原先顺序插入
则方案数可以分为两部分考虑:
2-1.环内排列，固定的环首不可排列，即((j-1)!)
2-2.剩下的(j-1)一个点位置未知，从未固定的(i+j-1)个点中选择
即(C(i+j-1,j-1))
所以就是(C(i+j-1,j-1)(j-1)!=frac{(i+j-1)!}{i!})
归纳一下，发现更形象的描述就是(egin{aligned}frac{n!}{prod_{i=1}^{m} (sum_{j=1}^i a_j)}end{aligned})
也就是每次除掉转移时的大小，将(n!)分成若干段，这似乎有利于理解下面的dp优化
3.计算(lcm)之积
考虑对于每个质因数计算其出现的幂次，注意这个幂次是对于(varphi)取模的
原先是求恰好包含(x^k)的方案数，得到的(dp)不好优化，考虑转换为:
求质因数(x)出现在答案里的幂次(ge k)的方案数(F_k)，答案就是(x^{sum F_k})
Solution 1
那么反向求解，令(dp_i)表示当前已经确定了(i)个点，没有出现(x^k)倍数大小的联通块
暴力转移，枚举(i)从所有(x^k ot|j)转移过来即可，单次求解复杂度为(O(n^2))，不可行
优化1:
考虑分解系数(frac{(i+j-1)!}{i!})，累前缀和，对于(j)为(x^k)倍数的情况枚举减掉
这样单次求解复杂度为(O(frac{n^2}{x^k}))，总复杂度为(O(n^2ln n))，且不好处理阶乘逆元
优化2:
不枚举(x^i)的倍数，直接再用一个前缀和数组(s_i)记录下来，让(s_i)从(s_{i-x^k})转移过来即可
如何将系数(frac{(i+j-1)!}{i!})分解？
每次(i+j)增大1，就多乘上一个(i+j)即可
当(s_{i})从(s_{i-x^k})转移过来时，需要补上(egin{aligned}prod_{j=i-x^k+1}^{i-1}jend{aligned})
也就是模拟了上面提到的把(n!)分段的过程
这样就去掉了阶乘逆元的求解
Tips:发现需要预处理(T_{i,j}=prod_{k=i}^j k)，可以滚动一下会快一点，内存为(O(n))
不同的(x^i)上限为(O(n))种，实际大概可能是(O(pi(n)log log n=frac{nlog log n}{log n}))?
因此复杂度为(O(n^2))
可以看到代码还是很简单的
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define Mod2(x) ((x<0)&&(x+=P2))
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)
const int N=7510;
int n,P,P2;
int mk[N];
int s[N],T[N],dp[N];
ll qpow(ll x,ll k=P-2) {
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
int main(){
scanf("%d%d",&n,&P),P2=P-1;
int ans=1;
rep(i,2,n) if(!mk[i]) {
for(int j=i+i;j<=n;j+=i) mk[j]=1;
for(int x=i;x<=n;x*=i) mk[x]=i;
}
rep(i,1,n) T[i]=1;
int r=1;
rep(i,1,n) r=1ll*r*i%P2;
rep(x,2,n) {
rep(j,1,n-x+1) T[j]=1ll*T[j]*(j+x-2)%P2;
// 滚动求解区间乘积
if(mk[x]<=1) continue;
dp[0]=1,s[0]=1;
int sum=1;
rep(i,1,n) {
s[i]=0;
if(i>=x) s[i]=1ll*s[i-x]*T[i-x+1]%P2;
dp[i]=sum-s[i]; Mod2(dp[i]);
s[i]=(1ll*s[i]*i+dp[i])%P2;
sum=(1ll*sum*i+dp[i])%P2;
}
ans=1ll*ans*qpow(mk[x],P2+r-dp[n])%P;
}
printf("%d
",ans);
}
Solution 2
为了便于表达，设满足条件为至少出现一个(x)的倍数
实际用min-max容斥确实比较好理解，设对于集合(S)，求其最大值
(egin{aligned} max lbrace S brace =sum_{Tsube S} (-1)^{|T|+1}minlbrace T braceend{aligned})
简要证明的话:
把(S)中的元素倒序排成一排分别为(S_i)
对于(S_1)即最大值，显然被计算一次
对于剩下的值(S_i(i>1))，则它作为最小值产生贡献意味着选的数都在(1-i)内，显然有(2^{i-1})次为奇数集合大小，(2^{i-1})为偶数集合大小，两部分抵消
要计算最大值为1的方案数，那么就要计算最小值为1的子集方案数
考虑强制一个子集中每一个环大小均为(x)的倍数，设选出了(i)个这样的环，总大小为(j)的方案数为(dp_{i,j})
则实际对答案的贡献还要考虑这样的子集出现的次数
考虑选择子集的位置，以及剩下的(n-j)个点任意排布，方案数应该为(C(n,j)cdot (n-j)!=frac{n!}{j!})
如果真的用(dp_{i,j})，复杂度显然太高，考虑(i)这一维的影响只在于系数((-1)^{|T|+1})，可以直接在转移过程中解决
因此可以直接记录大小(i)，从前面转移过来
(可以看到依然需要访问上面提到的(T_{i,j})，要滚动的话还会更难处理)
这样的(dp)状态有(frac{n}{x})种，转移为((frac{n}{x})^2)，最后统计复杂度为(O(frac{n}{x}))
实际上这样的复杂度已经足够了，是(O(sum frac{n^2}{i^2}approx n^2))
以下是(O(n^2))内存的代码
#include<bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for(int i=a;i<=b;++i)
#define drep(i,a,b) for(int i=a;i>=b;--i)
typedef long long ll;
const int N=7510;
int n,P,P2,T[N][N],mk[N],dp[N];
ll qpow(ll x,ll k) {
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
int Calc(int x){
int m=n/x,s=P2-1;
rep(i,1,m) {
s=1ll*s*T[(i-1)*x+1][i*x-1]%P2;
dp[i]=P2-s;
s=(1ll*s*i*x+dp[i])%P2;
}
int ans=0;
rep(i,1,m) ans=(ans+1ll*dp[i]*T[i*x+1][n])%P2;
return ans;
}
int main(){
scanf("%d%d",&n,&P),P2=P-1;
rep(i,2,n) if(!mk[i]) {
for(int j=i;j<=n;j+=i) mk[j]=1;
for(int j=i;j<=n;j*=i) mk[j]=i;
}
rep(i,1,n+1){
T[i][i-1]=1;
rep(j,i,n) T[i][j]=1ll*T[i][j-1]*j%P2;
}
int ans=1;
rep(i,2,n) if(mk[i]>1) ans=ans*qpow(mk[i],Calc(i))%P;
printf("%d
",ans);
}
最终优化
可以看到，这个做法和Sol1的转移有十分的相同之处，因此考虑用同样的方法优化掉
但是预处理系数的部分依然是(O(n^2))的，如何解决呢?
1.线段树大法
预处理复杂度为(O(n))，查询复杂度为(O(log n))，总复杂度(O(nlog^2 n))
空间复杂度为(O(n))
Oh这个代码是ZKW线段树
#include<bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for(int i=a;i<=b;++i)
typedef long long ll;
const int N=7510;
int n,P,P2;
struct Tree {
int s[N<<2],bit;
void Build() {
for(bit=1;bit<=n;bit<<=1);
rep(i,1,n) s[i+bit]=i;
for(int i=bit;i>=1;--i) s[i]=1ll*s[i<<1]*s[i<<1|1]%P2;
}
int Que(int l,int r){
if(l>r) return 1;
if(l==r) return l;
int res=1;
for(l+=bit-1,r+=bit+1;l^r^1;l>>=1,r>>=1){
if(~l&1) res=1ll*res*s[l^1]%P2;
if(r&1) res=1ll*res*s[r^1]%P2;
}
return res;
}
} T;
int mk[N];
ll qpow(ll x,ll k) {
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
int dp[N];
int Calc(int x){
int m=n/x,s=P2-1;
rep(i,1,m) {
s=1ll*s*T.Que((i-1)*x+1,i*x-1)%P2;
dp[i]=P2-s;
s=(1ll*s*i*x+dp[i])%P2;
}
int ans=0;
rep(i,1,m) ans=(ans+1ll*dp[i]*T.Que(i*x+1,n))%P2;
return ans;
}
int main(){
scanf("%d%d",&n,&P),P2=P-1;
T.Build();
rep(i,2,n) if(!mk[i]) {
for(int j=i;j<=n;j+=i) mk[j]=1;
for(int j=i;j<=n;j*=i) mk[j]=i;
}
int ans=1;
rep(i,2,n) if(mk[i]>1) ans=ans*qpow(mk[i],Calc(i))%P;
printf("%d
",ans);
}
2.模数分解+CRT(Chinese Remainder Theory)
分解后，可以用模逆元处理，然后就直接做，最后CRT合并一下，其实我并不会实现。。。
3.猫树(嘿嘿)
这是一个(O(nlog n))预处理，(O(1))查询的数据结构，空间复杂度为(O(nlog n))
因此复杂度为(O(nlog n))
#include<bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for(int i=a;i<=b;++i)
#define drep(i,a,b) for(int i=a;i>=b;--i)
typedef long long ll;
const int N=7510;
int n,P,P2;
struct SuckCat{
int s[14][8200],Log[8200],bit;
void Build() {
for(bit=1;bit<=n;bit<<=1);
rep(i,2,bit) Log[i]=Log[i>>1]+1;
for(int l=1,d=0;l*2<=bit;l<<=1,d++){
for(int i=l;i<=bit;i+=l*2){
s[d][i-1]=i-1,s[d][i]=i;
drep(j,i-2,i-l) s[d][j]=1ll*s[d][j+1]*j%P2;
rep(j,i+1,i+l-1) s[d][j]=1ll*s[d][j-1]*j%P2;
}
}
}
int Que(int l,int r){
if(l>r) return 1;
if(l==r) return l;
int d=Log[l^r];
return 1ll*s[d][l]*s[d][r]%P2;
}
} T;
int mk[N];
ll qpow(ll x,ll k) {
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
int dp[N];
int Calc(int x){
int m=n/x,s=P2-1;
rep(i,1,m) {
s=1ll*s*T.Que((i-1)*x+1,i*x-1)%P2;
dp[i]=P2-s;
s=(1ll*s*i*x+dp[i])%P2;
}
int ans=0;
rep(i,1,m) ans=(ans+1ll*dp[i]*T.Que(i*x+1,n))%P2;
return ans;
}
int main(){
scanf("%d%d",&n,&P),P2=P-1;
T.Build();
rep(i,2,n) if(!mk[i]) {
for(int j=i;j<=n;j+=i) mk[j]=1;
for(int j=i;j<=n;j*=i) mk[j]=i;
}
int ans=1;
rep(i,2,n) if(mk[i]>1) ans=ans*qpow(mk[i],Calc(i))%P;
printf("%d
",ans);
}