自动搬运

查看原文

来自洛谷，原作者为

	NaCly_Fish 北海虽赊，扶摇可接。

搬运于2025-08-24 23:09:52，当前版本为作者最后更新于2025-04-23 16:09:48，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

零

尝试用尽可能少的前置知识来推导。
如果你对这种大力推导的方法比较感兴趣，那我推荐你来看看。

看到这种问题，一般的想法是对于 $n$ 种素数（也是每个维度）独立处理，然后把答案再合并起来。

可是这里没那么简单，比如对于子问题 1，我们可以用 $n+1$ 元生成函数写出答案的最朴素的表达式：

$$\left[y^k \prod_{i=1}^n x_i^{a_i} \right] \left(\prod_{\bold i}\left( 1+y \prod_{j=1}^n x_j^{i_j}\right)\right)\left(\prod_{\bold i \leq \bold b}\left( 1+y\prod_{j=1}^n x_j^{i_j}\right) \right)^{-1} $$

其中 $\bold i = (i_1,\cdots,i_n)$，$\bold b=(b_1,\cdots,b_n)$。式子中用 $y$ 计量使用的正整数数量，$x_1,\cdots,x_n$ 分别计量 $n$ 种素数的幂次。

因为这个对 $y$ 的限制很麻烦，如果能尝试把它干掉，那么对于 $x_1,\cdots,x_n$ 大概就能独立处理了吧？

---

一

不过原问题还是太复杂，可以先考虑这样一个特别简化后的情况：不必满足 $M$ 的限制，且 $n=1$ 时的子问题 1。

此时要求的就是「将一个正整数拆分为若干互异自然数之和」的方案数，我们容易写出答案的生成函数为：

$$[y^3]\prod_{i=0}^\infty (1+y x^i)=\frac 16\frac{1}{(1-x)^3}-\frac 12\frac{1}{(1-x)(1-x^2)}+\frac 13\frac{1}{(1-x^3)} $$

等号左边就是我们刚才写过的，但右边有什么用，又是怎么得到的呢？

这个东西厉害的地方在于它可以很容易扩展至多元的情况，即：

$$[y^3] \prod_{i_1,\cdots,i_n}\left(1+y\prod_{j=1}^n x_j^{i_j} \right)=\frac 16\prod_{j=1}^n \frac{1}{(1-x_j)^3}-\frac 12\prod_{j=1}^n \frac{1}{(1-x_j)(1-x_j^2)}+\frac13\prod_{j=1}^n \frac{1}{(1-x_j^3)} $$

这样就可以轻松解决不限制 $n$ 时 $k=3$ 的情况了（但还是不要求满足 $M$ 的限制）。经过这样一通操作，我们确实将这 $n$ 个维度独立开来，可以分别计算了。再提取其 $x_1^{a_1} \cdots x_n^{a_n}$ 系数即为：

$$\frac 16\prod_{j=1}^n\binom{a_j+2}{2}-\frac 12\prod_{j=1}^n (\lfloor a_j/2\rfloor +1)+\frac 13 \prod_{j=1}^n[3\mid a_j] $$

非常的简单方便。

对于刚才的第二个问题，答案也很简单：求导即可。

设

$$F(y;x_1,\cdots,x_n) = \prod_{i_1,\cdots,i_n}\left(1+y \prod_{j=1}^nx_j^{i_j}\right) $$

那么 $[y^k]F(y) = F^{(k)}(0)/k!$（Taylor 展开）。用经典的对数求导法即得

$$F'(y)=F(y)\sum_{i_1,\cdots,i_n} \frac{\prod_{j=1}^n x_j^{i_j}}{1+y \prod_{j=1}^n x_j^{i_j}} $$

再对等式两边同求 $(k-1)$ 阶导（使用 Leibniz 公式），然后代入 $y=0$：

$$F^{(k)}(0)=\sum_{i=0}^{k-1}\binom{k-1}{i} \frac{(-1)^i i!}{\prod_{j=1}^n(1-x_j^{i+1})} F^{(k-i-1)}(0) $$

边界条件自然是 $F(0)=1$，根据此式递推计算即可。

对于子问题 2 当然也很简单，把 $(-1)^i$ 去掉就好了。

---

二

如果要对子问题 1 引入不为 $M$ 约数的限制，也是简单的。设

$$\hat F(y) \prod_{\bold i \leq \bold b}\left(1+y \prod_{j=1}^n x_j^{i_j} \right)= \prod_{i_1,\cdots,i_n}\left(1+y \prod_{j=1}^nx_j^{i_j}\right) $$

按上面的方法算一遍，同时我们简记

$$U_i=(-1)^{i+1}\prod_{j=1}^n \frac{1}{1-x_j^i} \ , \ V_i=U_i \prod_{j=1}^n(1-x_j^{i(b_j+1)}) $$

那么就有

$$\hat F^{(k)}(0) = \sum_{i=0}^{k-1}\binom{k-1}{i}i! (U_{i+1}-V_{i+1})\hat F^{(k-i-1)}(0) $$

再设 $W_i=U_i-V_i$。

很显然 $\hat F^{(k)}(0)$ 是由 $W_1,\cdots,W_k$ 构成的多项式表示的。进一步观察发现，其中的每一项 $W$ 的下标之和（有重复的也要再记）都等于 $k$。

比如

$$\hat F^{(4)}(0)=W_1^4+6W_1^2W_2+3W_2^2+8W_1W_3+6W_4 $$

那这显然就和分拆数问题对应起来了，$\hat F^{(k)}(0)$ 中仅有 $P_k$ 项（$P_k$ 即为将 $k$ 拆分为若干无序正整数之和的方案数）。在下一部分中，我们将给出每一项的系数。

有同学看到这里可能按捺不住了：「Burnside 引理不就直接秒了吗？」但是 Burnside 引理还是太高级了，我们要用更朴素的方法。

---

三

我们把 $\hat F^{(k)}(0)$ 的递归式再整理一下：

$$\frac{ \hat F^{(k)}(0)}{k!}=\frac 1k\sum_{i=0}^{k-1}W_{i+1} \frac{F^{(k-i-1)}(0)}{(k-i-1)!} $$

设 $G_k=\hat F^{(k)}(0)/k!$，并考虑 $G_k$ 中的某一项 $W_1,\cdots,W_k$ 分别出现了 $c_1,\cdots,c_k$ 次。

只关注其中一项，比如 $k=4$ 时 $c=(2,1,0,0)$ 的这一项系数是 $1/4$，它是这么来的：

$$\frac{1}{4}W_1\frac{1}{3}W_1\frac 12W_2+\frac14W_1\frac 13W_2\frac{1}{1}W_1+\frac 14W_2\frac 12W_1\frac11W_1=\frac 14W_1^2W_2 $$

这里就是枚举了三种 $W_i$ 在递归中乘的顺序，然后对应乘上和式前面 $1/k$ 的系数。

扩展到一般的情况，这一项的系数就是对于全部 $\binom{c_1+\cdots +c_k}{c_1,\cdots,c_k}$ 种排列，求其「前缀和序列之积的倒数」之和。
再仔细想想，可以将其对应到「确定 $k$ 阶排列中置换环结构后的方案计数」上去（暂时没有很好的代数推导方法，待补充）。

所以这个式子实际算的就是：

$$\prod_{i=1}^k \frac{1}{c_i ! i^{c_i}}$$

现在有了每一项的系数，就可以直接展开 $W_i$ 计算了。

---

四

将所有 $(U_i - V_i)$ 的乘积全部展开来，每一项只观察 $n$ 维中的一部分，它是这样的形式（略去了其整体前面的系数）：

$$\prod_{i=1}^k \frac{(1-z^{i(b+1)})^{d_i}}{(1-z^i)^{c_i}} $$

其中 $d_i \leq c_i$，$b$ 就是这一维中对应的 $M$ 的限制。

将分子展开，最终就只需要多次给出 $N$ 并查询：

$$[z^N] \prod_{i=1}^k \frac{1}{(1-z^i)^{c_i}}$$

设 $L= \text{lcm}\{ i^{[c_i > 0]}\}$，那么上式可以表示为 $P_{N\bmod L}(\lfloor N/L \rfloor)$，其中 $P_i(x)$ 是一组 $k$ 次多项式。

那么只需要计算出 $\prod_i (1-z^i)^{-c_i}$ 的前 $L(k+1)$ 项系数，就可以使用线性插值求单点的方法，快速预处理出所有将要用到的点值。

---

五

参考代码如下，加了一些注释，可以补充文中没描述清楚的地方：

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#define N 203
#define p 1000000021
#define ll long long
using namespace std;

inline int power(int a,int t){
	int res = 1;
	while(t){
		if(t&1) res = (ll)res*a%p;
		a = (ll)a*a%p;
		t >>= 1;
	}
	return res;
}

int fac[N],ifac[N],inv[N];

void init(int n){
	fac[0] = ifac[0] = inv[1] = 1;
	for(int i=1;i<=n;++i) fac[i] = (ll)fac[i-1]*i%p;
	ifac[n] = power(fac[n],p-2);
	for(int i=n-1;i;--i) ifac[i] = (ll)ifac[i+1]*(i+1)%p;
	for(int i=2;i<=n;++i) inv[i] = (ll)fac[i-1]*ifac[i]%p;
}

inline int binom(int n,int m){
	return (ll)fac[n]*ifac[m]%p*ifac[n-m]%p;
}

int gcd(int a,int b){ return b==0?a:gcd(b,a%b); }
int lcm(int a,int b){ return a*b/gcd(a,b); }

int interpolate(const int *f,int n,int x){ // give n-degree poly f values f(0) , ... , f(n), return f(x).
	static int pre[N],suf[N];
	n += 1,x += 1;
	pre[0] = suf[n+1] = 1;
	for(int i=1;i<=n;++i) pre[i] = (ll)pre[i-1]*(x-i)%p;
	for(int i=n;i;--i) suf[i] = (ll)suf[i+1]*(x-i)%p;
	int ans = 0,tmp;
	for(int i=1;i<=n;++i){
		tmp = (ll)f[i-1]*ifac[i-1]%p*ifac[n-i]%p;
		if((n-i)&1) tmp = -tmp;
		ans = (ans + (ll)tmp*pre[i-1]%p*suf[i+1])%p;
	}
	return (ans+p)%p;
}

ll a[N],b[N];
int c[N],tc[N],val[53][N];
int n,k,ans1,ans2;
int delt;

void dfs(int dep,int cf){ // expand \prod_{i=1}^k (U_i - V_i)^{c_i}
	if(dep > k){
		int flag = 0;
		for(int i=1;i<=k;++i){
			cf = (ll)cf*binom(c[i],tc[i])%p;
			flag += tc[i];
		}
		if(flag&1) cf = -cf;
		int f[101]; // the numerator. f(z) = \prod_{i=1}^k (1-z^i)^{tc_i}
		memset(f,0,sizeof(f));
		f[0] = 1;
		int deg = 0;
		for(int i=1;i<=k;++i)
		for(int t=1;t<=tc[i];++t){
			for(int j=deg+i;j>=i;--j)
				f[j] = (f[j] - f[j-i])%p;
			deg += i;
		}
		int tmp = 1;
		for(int i=1;i<=n;++i){
			int sum = 0;
			for(int j=0;j<=deg;++j)
				sum = (sum + (ll)f[j]*val[i][j])%p;
			tmp = (ll)tmp*sum%p;
		}
		flag = 0;
		for(int i=1;i<=k;++i) flag += (i+1)*c[i];
		ans2 = (ans2 + (ll)tmp*cf)%p;
		if(flag&1) tmp = p-tmp;
		ans1 = (ans1 + (ll)tmp*cf)%p;
		return;
	}
	for(int i=0;i<=c[dep];++i){
		tc[dep] = i;
		dfs(dep+1,cf);
	}
}

void solve(){
	int d = 1,cf = 1,m = 0;
	for(int i=1;i<=k;++i){
		if(c[i] > 0) d = lcm(d,i);
		m += c[i];
		cf = (ll)cf*ifac[c[i]]%p*power(inv[i],c[i])%p;
	}
	static int f[10003],g[2000][N];
	memset(f,0,sizeof(f));
	memset(val,0,sizeof(val));
	int L = d*(m+1);
	f[0] = 1;
	for(int i=1;i<=k;++i)
	for(int t=1;t<=c[i];++t)
	for(int j=i;j<=L;++j)
		f[j] = (f[j] + f[j-i])%p;

	for(int i=0;i<d;++i)
	for(int j=0;j<=m;++j)
		g[i][j] = f[j*d+i];

	for(int i=1;i<=n;++i)
	for(int j=0;j<=k;++j){
		if(a[i] < j*(b[i]+1)) break;
		ll x = a[i] - j*(b[i]+1);
		val[i][j] = interpolate(g[x%d],m,(x/d)%p);
        // val[i][j] = [z^{a_i - j(b_i+1)}] \prod_{t=1}^k 1/(1-z^t)^{c_t}
    }
	dfs(1,cf);
}

void partition(int sum,int dep){ // find all partition of k
	if(sum == k){
		solve();
		return;
	}
	if(dep > k-sum) return;
	for(int i=0;sum+i*dep<=k;++i){
		c[dep] = i;
		partition(sum+i*dep,dep+1);
	}
	c[dep] = 0;
}


int main(){
	scanf("%d%d",&n,&k);
	for(int i=1;i<=n;++i) scanf("%lld",&a[i]);
	for(int i=1;i<=n;++i) scanf("%lld",&b[i]);
	init(103);
	partition(0,1);
	printf("%d %d",(ans1+p)%p,(ans2+p)%p);
	return 0;	
}