P11739 朴素推导

· · 题解

零

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

看到这种问题，一般的想法是对于 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;   
}