集训队互测 2015 普罗达科特
令 \(N=\prod p_i^{a_i},M=\prod p_i^{b_i}\),\(p\) 为两两不同的素数,\(1\le i\le n\)。求有多少本质不同的大小为 $m $ 的不可重集1和可重集2 \(S\) 使得 \(S\) 的元素乘积为 \(N\) 且每个元素都不整除 \(M\)。\(m\le 25,n\le 50\).
先讨论不可重集的问题,然后先不管这个 \(M\)。
我们看看怎么把问题给化简。我们对相等关系进行容斥。假如我们暴力枚举每两个元素之间是否钦定相等关系,那么复杂度是很不对的。但是我们实际上只关心这些等价类的大小,所以我们可以枚举划分。而一个划分 \(\pi:p_1+\dots+p_k=n\) 的权值是所有能用这个划分表示的图的 \(-1\) 的边数次方之和。我们设 \(f_x\) 表示大小为 \(x\) 的连通图的 \(-1\) 的边数次方之和,于是我们可以容斥 DP,枚举 \(1\) 所在连通块大小。但是你发现由于在点数 \(>1\) 时任意图的权值之和 \(\sum_i (-1)^i\binom{x(x-1)/2}{i}=0\),没有贡献,所以我们实际上可以化简成为 \(f_x=-(n-1)f_{x-1}\),于是可以得到 \(f_x=(-x+1)!=(-1)^{x-1}(x-1)!\)。所以对于这样的一个划分,设每个元素出现了 \(q_i\) 次,那么可以得到其系数(化简后)为
然后化简可重集的问题。我们发现可重集的问题是什么呢?相当于就是在置换群作用下的本质不同方案个数!我们发现不动点相当于每个轮换都形成一个等价类。于是我们枚举轮换大小划分,其中每个划分方案的系数为
我们发现好想两个 \(g\) 好像就差一个 \((-1)^{m-k}\)。有趣的。
于是现在两个问题都转化成了同一个问题:就是给定一些等价类,给每个等价类弄组权值,使得满足条件。
现在我们开始要解决如下的问题:给定了一个整数拆分(即每个等价类大小),将带标号元素划分到集合中,然后满足题目限制。
我们还是很讨厌这个 \(M\) 的限制。\(M\) 的限制具体就是说,设每个等价类在这些素数上取的值分别为 \(d_1,\dots,d_n\),那么至少要有一个 \(d>b\)。我们考虑容斥。由于大小相同的等价类本质相同,所以我们枚举对于大小为 \(x\) 的等价类(有 $ q_i$ 个)中有 \(c_i\) 个被钦定不满足限制,于是得到容斥系数 \(\prod_i \binom{q_i}{c_i}(-1)^{c_i}\)。你发现这个 \(c\) 也确实是可以暴力枚举的,需要枚举的 \(q,c\) 组合的数量级 \(Q_m\) 大概在 \(10^5\)。
现在我们就成功把素因子给独立出来了。对于每一个素因子,我们枚举每一个等价类大小 \(j\),我们需要满足其中有 \(c_j\) 个的取值在 \(b_i\) 之内。考虑使用生成函数。我们对于一个没有被钦定的等价类,它的生成函数显然为 \(\frac{1}{1-z^j}\)。而对于一个被钦定的等价类,我们补集转换,减去大于的方案数。所以生成函数为 \(\frac{1-z^{j(b_i+1)}}{1-z^j}\)。所以可以得到总生成函数数为
如果摁做复杂度会爆。对于这些分子,注意到这些 \(z\) 都是 \(b_i+1\) 的倍数,所以它本质是一个关于 \(x^{b_i+1}\) 的 \(\sum jc_j\le m\) 次多项式 \(T\),所以暴力展开的复杂度是 \(O(m^2)\)。对于分母看上去不大好做,但是我们实际上只需要对于一些 \(l\) 求出 \([z^{a_i-l(b_i+1)}]\prod \frac{1}{1-p_j}\)。有一个非常劲爆的引理:
\(\sum_{i=1}^{m} a_ix_i=N\) 的非负整数解的组数 \(F_N=[z^N]\prod_i \frac{1}{1-z^{a_i}}\),记 \(L=\operatorname{lcm}(a_i)\),那么 \(F_{kL+r}\) 是关于 \(k\) 的 \(m\) 次多项式。
于是我们就能处理这个分母 \(S\) 了。我们 \(O(P_mnm^2+S_mm^2+Q_mnm))\) 求出所有我们需要的值。\(S=\sum L\)。
总复杂度 \(O(P_mnm^2+S_mm^2+Q_m(m^2+nm))\)。可能复杂度写的有些细节问题。
#include<bits/stdc++.h>
#define int long long
#define rep(i,a,b) for(int i=(a);i<=(b);i++)
#define per(i,a,b) for(int i=(a);i>=(b);i--)
#define fi first
#define se second
#define eb emplace_back
#define popc __builtin_popcount
#define sgn(x) ((x)&1?-1:1)
using namespace std;
typedef long long ll;
typedef pair<int,int> pii;
typedef vector<int> vi;
typedef vector<pii> vp;
typedef unsigned long long ull;
typedef long double ld;
int read() {
int x=0,w=1; char c=getchar();
while(!isdigit(c)) {if(c=='-') w=-1; c=getchar();}
while(isdigit(c)) {x=x*10+c-'0'; c=getchar();}
return x*w;
}
namespace Sol {
const int N=50009,mod=1e9+21;
int inv[N],fac[N],ifac[N];
int ksm(int x,int y,int r=1) {
for(;y;y>>=1,x=x*x%mod) if(y&1) r=r*x%mod;
return r;
}
void pre(int n) {
inv[0]=inv[1]=fac[0]=ifac[0]=1;
rep(i,1,n) fac[i]=fac[i-1]*i%mod;
ifac[n]=ksm(fac[n],mod-2);
per(i,n-1,1) ifac[i]=ifac[i+1]*(i+1)%mod;
rep(i,2,n) inv[i]=ifac[i]*fac[i-1]%mod;
}
struct Poly {
int deg; vi f;
Poly(int x=0) {deg=x; f.resize(deg+1);}
Poly(vi p) {f=p; deg=(int)p.size()-1;}
int& operator [] (int x) {return f[x];}
int itp(int x,int y=0) { //当表示点值时的 lagrange interpolation
if(x<=deg) return f[x]; static int s[59], t[59];
s[0]=1; rep(i,1,deg) s[i]=s[i-1]*(x-i+1)%mod;
t[deg]=1; per(i,deg-1,0) t[i]=t[i+1]*(x-i-1)%mod;
rep(i,0,deg) y=(y+sgn(deg-i)*s[i]%mod*t[i]%mod*ifac[i]%mod*ifac[deg-i]%mod*f[i])%mod;
return (y+mod)%mod;
}
};
Poly operator * (Poly &a,Poly &b) {
Poly c(a.deg+b.deg);
rep(i,0,a.deg) rep(j,0,b.deg) c[i+j]=(c[i+j]+a[i]*b[j])%mod;
return c;
}
Poly invs(Poly a,int deg) {
Poly b(deg);
b[0]=ksm(a[0],mod-2);
rep(i,1,deg) {
rep(j,1,min(a.deg,i)) b[i]=(b[i]+a[j]*b[i-j])%mod;
b[i]=mod-b[i]*b[0]%mod; if(b[i]>=mod) b[i]-=mod;
} return b;
}
int n,m,a[N],b[N],l,q[N],c[N],g0,g1,s[59][59],ans0,ans1;
Poly S,T,h[N];
void cans() {
T=Poly(0); T[0]=1; int res=1;
rep(x,1,m) if(c[x]) {
Poly r(x); r[0]=1, r[x]=mod-1; int y=c[x];
for(;y;y>>=1,r=r*r) if(y&1) T=T*r;
}
rep(i,1,n) {
int sum=0;
rep(j,0,T.deg) sum=(sum+T[j]*s[i][j])%mod;
res=res*sum%mod;
}
rep(x,1,m) res=res*ifac[c[x]]%mod*ifac[q[x]-c[x]]%mod*sgn(c[x]);
res=(res+mod)%mod;
ans0=(ans0+res*g0)%mod, ans1=(ans1+res*g1)%mod;
}
void dfs_exc(int x) {
if(x==m+1) {cans(); return;}
rep(i,0,q[x]) c[x]=i, dfs_exc(x+1);
}
void calc() {
S=Poly(0); S[0]=1; l=1;
rep(x,1,m) if(q[x]) l=l*x/__gcd(l,x);
rep(x,1,m) if(q[x]) {
Poly r(x); r[0]=1, r[x]=mod-1; int y=q[x];
for(;y;y>>=1,r=r*r) if(y&1) S=S*r;
}
S=invs(S,l*(m+1)-1);
rep(i,1,n) rep(x,0,m) if(a[i]/(b[i]+1)>=x) {
int y=a[i]-x*(b[i]+1), k=y/l, r=y-k*l;
Poly g(m); rep(j,0,m) g[j]=S[j*l+r];
s[i][x]=g.itp(k%mod);
}
g0=g1=1;
rep(x,1,m) rep(j,1,q[x]) {
g0=(g0*fac[x-1]%mod*ifac[x]%mod*sgn(x-1)+mod)%mod;
g1=g1*fac[x-1]%mod*ifac[x]%mod;
}
dfs_exc(1);
}
void dfs_div(int x,int y) {
if(x==m+1) {if(y==m) calc(); return;}
rep(i,0,(m-y)/x) q[x]=i, dfs_div(x+1,y+x*i);
}
void main() {
n=read(), m=read(), pre(50000);
mt19937 rnd(time(0));
rep(i,1,n) a[i]=read();
rep(i,1,n) b[i]=read();
dfs_div(1,0);
printf("%lld %lld\n",ans0,ans1);
}
}