集训队互测 2015 普罗达科特-526互联

令 $N=\prod p_i^{a_i},M=\prod p_i^{b_i}$，$p$ 为两两不同的素数，$1\le i\le n$。求有多少本质不同的大小为 $m $ 的不可重集¹和可重集² $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=\frac{(-1)^{m-k}}{\prod q_i!\prod p_i} \]

然后化简可重集的问题。我们发现可重集的问题是什么呢？相当于就是在置换群作用下的本质不同方案个数！我们发现不动点相当于每个轮换都形成一个等价类。于是我们枚举轮换大小划分，其中每个划分方案的系数为

\[g=\frac{1}{\prod q_i!\prod p_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}$。所以可以得到总生成函数数为

\[\prod_{j} \frac{(1-z^{j(b_i+1)})^{c_j}}{(1-z^j)^{q_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);
}

}