BZOJ4126 国王奇遇记加强版之再加强版 <多项式插值>

Problem

国王奇遇记加强版之再加强版

$\mathrm{Time\;Limit:\;15\;Sec}$
$\mathrm{Memory\;Limit:\;128\;MB}$

Description

Input

共一行，包括两个正整数$N$和$M$。

Output

共一行，为所求表达式的值对$10^9+7$取模的值。

Sample Input

5 3

Sample Output

36363

Hint

$1\le N\le10^9$，$1\le m\le5\times10^5$

Source

$\mathrm{By\;wjy1998}$

标签：多项式插值

Solution

好题，看了$\mathrm{Miskcoo's\;Space}$中的题解才懂。以下题解全部部分摘自特殊多项式在整点上的线性插值方法和BZOJ-3157. 国王奇遇记。

1. 多项式整点插值

观察二项式系数$\binom{x}{m}=\frac{x(x-1)(x-2)\cdots(x-m+1)}{m!}$，其为一个$x$的$m$次多项式。对于$\binom{x}{0},\binom{x}{1},\cdots,\binom{x}{m}$，由于其次数互不相同，故其线性无关。可以发现这$m+1$个多项式是$m$次多项式线性空间$P_m[x]$的一组基。
于是对于$\forall F(x)\in P_m[x]$，其一定可以表示为这$m+1$个多项式的线性组合，即

$$F(x)=\sum_{i=0}^{m}\binom{x}{i}a_i\tag{1}$$

由于$i>x$时，$\binom{x}{i}=0$，于是当$x>m$时，有

$$F(x)=\sum_{i=0}^{x}\binom{x}{i}a_i\tag{2}$$

根据二项式定理，可知$\sum_{i=0}^{m}(-1)^i\binom{m}{i}=[m=0]$，应用其对$(2)$进行二项式反演，得

$$a_i=\sum_{j=0}^{i}(-1)^{i-j}\binom{i}{j}F(j)\tag{3}$$

将$(3)$带入$(1)$，得

$$\begin{aligned} F(x) &=\sum_{i=0}^{m}\binom{x}{i}\sum_{j=0}^{i}(-1)^{i-j}\binom{i}{j}F(j)\\ &=\sum_{j=0}^{m}F(j)\sum_{i=j}^{m}\binom{x}{i}(-1)^{i-j}\binom{i}{j}\\ &=\sum_{j=0}^{m}F(j)\sum_{i=0}^{m-j}(-1)^{i}\binom{x}{i+j}\binom{i+j}{j}\\ \end{aligned}$$

化简一下二项式系数

$$\begin{aligned} \binom{x}{i+j}\binom{i+j}{j} &=\frac{x!}{(i+j)!(x-i-j)!}\times\frac{(i+j)!}{i!j!}\\ &=\frac{x!}{i!j!(x-i-j)!}\\ &=\frac{x!}{j!(x-j)!}\times\frac{(x-j)!}{i!(x-i-j)!}\\ &=\binom{x}{j}\binom{x-j}{i}\\ \end{aligned}$$

于是

$$\begin{aligned} F(x) &=\sum_{j=0}^{m}F(j)\sum_{i=0}^{m-j}(-1)^{i}\binom{x}{i+j}\binom{i+j}{j}\\ &=\sum_{j=0}^{m}F(j)\sum_{i=0}^{m-j}(-1)^{i}\binom{x}{j}\binom{x-j}{i}\\ &=\sum_{j=0}^{m}\binom{x}{j}F(j)\sum_{i=0}^{m-j}(-1)^{i}\binom{x-j}{i}\\ \end{aligned}$$

即

$$F(x)=\sum_{j=0}^{m}\binom{x}{j}F(j)\sum_{i=0}^{m-j}(-1)^{i}\binom{x-j}{i}\tag{4}$$

将后面的部分进一步化简，即化简形如$\sum_{i=0}^{n}(-1)^i\binom{x}{i}$的式子。
根据上指标反转$\binom{p}{q}=(-1)^q\binom{q-p-1}{q}$，可将其化简

$$\begin{aligned} \sum_{i=0}^{n}(-1)^i\binom{x}{i} &=\sum_{i=0}^{n}\binom{i-x-1}{i}\\ &=\binom{-x-1+0}{0}+\binom{-x-1+1}{1}+\cdots+\binom{-x-1+n}{n}\\ &=\binom{-x-1+1}{0}+\binom{-x-1+1}{1}+\cdots+\binom{-x-1+n}{n}\\ &=\binom{-x-1+2}{1}+\binom{-x-1+2}{2}+\cdots+\binom{-x-1+n}{n}\\ &=\binom{-x-1+n+1}{n}\\ &=\binom{n-x}{n}\\ \end{aligned}$$

带入$(4)$得

$$\begin{aligned} F(x) &=\sum_{j=0}^{m}\binom{x}{j}F(j)\sum_{i=0}^{m-j}(-1)^{i}\binom{x-j}{i}\\ &=\sum_{j=0}^{m}\binom{x}{j}\binom{m-x}{m-j}F(j)\\ \end{aligned}$$

将后面的$\binom{x-j-1}{m-x}$再次上指标反转得

$$F(x)=\sum_{j=0}^{m}(-1)^{m-j}\binom{x}{j}\binom{x-j-1}{m-j}F(j)$$

即

$$F(x)=\sum_{i=0}^{m}(-1)^{m-i}\binom{x}{i}\binom{x-i-1}{m-i}F(i)$$

故对于$m$次多项式$F(x)$，当$x>m$时，如果能得到$F(0)\sim F(m)$的值，可以求出$F(x)$的值。
具体地，两个二项式系数的乘积为

$$\begin{aligned} \binom{x}{i}\binom{x-i-1}{m-i} &=\frac{x!}{i!(x-i)!}\times\frac{(x-i-1)!}{(m-i)!(x-m-1)!}\\ &=\frac{x!}{i!(x-i)(m-i)!(x-m-1)!}\\ &=\frac{x(x-1)(x-2)\cdots(x-m)}{x-i}\times\frac{i!}{(m-i)!}\\ &=x(x-1)\cdots(x-i+1)\times(x-i-1)(x-i-2)\cdots(x-m)\times\frac{i!}{(m-i)!}\\ &=\frac{i!}{(m-i)!}\prod_{p=x-m}^{x-i-1}p\prod_{q=x-i+1}^{x}q\\ \end{aligned}$$

对于分式部分，可以预处理阶乘及其逆元以$O(1)$计算，对于两个乘式，可以预处理$(x-m)\sim x$的前缀积和后缀积以$O(1)$计算。总时间复杂度为$O(m)$。

2. 幂和

记$S_m(n)=\sum_{i=1}^{n}i^mm^i$，求出$m$较小时的通项，瞎猜发现当$m>1$时，$S_m(n)$一定有如下形式：

$$S_m(n)=m^nF_m(n)-F_m(0)$$

其中$F_m(n)$是一个$m$次多项式。根据前面多项式整点线性插值，只需求出$F_m(0),F_m(1),\cdots,F_m(m)$即可$O(m)$求得$F_m(n)$。
注意到在$S_m(i)-S_m(i-1)$中$F_m(0)$被消去，于是可以得到$F_m$的递推式：

$$S_m(i)-S_m(i-1)=m^iF_m(i)-m^{i-1}F_m(i-1)=i^mm^i\\ F_m(i)=\frac{F_m(i-1)}{m}+i^m$$

于是可以将$F_m(i)$表示为$k_iF_m(0)+b_i$的形式，然而还无法求出$F_m(0)$。
根据多项式插值时推导的公式，当$x>m$时，有

$$F_m(x)=\sum_{i=0}^{m}(-1)^{m-i}\binom{x}{i}\binom{x-i-1}{m-i}F_m(i)$$

将$m+1$作为$x$带入得

$$\begin{aligned} F_m(m+1) &=\sum_{i=0}^{m}(-1)^{m-i}\binom{m+1}{i}\binom{m-i}{m-i}F_m(i)\\ &=\sum_{i=0}^{m}(-1)^{m-i}\binom{m+1}{i}F_m(i)\\ \end{aligned}$$

将$F(m+1)$移到右边得

$$\sum_{i=0}^{m}(-1)^{m-i}\binom{m+1}{i}F_m(i)=0$$

将$\forall i\in[1,n+1],F_m(i)=k_iF_m(0)+b_i$带入即可解出$F_m(0)$，复杂度$O(m)$。然后通过前面$O(m)$插值可求出$F_m(n)$，带入$S_m(n)=m^nF_m(n)-F_m(0)$即可求得$S_m(n)$。总时间复杂度$O(m)$。
注意由于$m=1$时不满足性质，需要单独计算；此外还需暴力计算$n\le m$的情况。

简单版之再简单版见BZOJ3157，简单版见BZOJ3516。

Code

#include <bits/stdc++.h>
#define MAX_M 500000
#define P 1000000007
using namespace std;
typedef long long lnt;
template <class T> inline void read(T &x) {
    x = 0; int c = getchar(), f = 1;
    for (; !isdigit(c); c = getchar()) if (c == 45) f = -1;
    for (; isdigit(c); c = getchar()) (x *= 10) += f*(c-'0');
}
int n, m; lnt F[MAX_M+5];
lnt pw[MAX_M+5], fac[MAX_M+5], inv[MAX_M+5];
bool NotPri[MAX_M+5]; vector <int> pri;
lnt Pow(lnt x, lnt k) {
    lnt ret = 1;
    for (; k; k >>= 1, x = x*x%P)
        if (k&1) ret = ret*x%P;
    return ret;
}
void init(int n) {
    fac[0] = inv[0] = inv[1] = pw[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!NotPri[i]) pri.push_back(i), pw[i] = Pow(i, m);
        for (int j = 0; j < (int)pri.size(); j++) {
            if (i*pri[j] > n) break;
            NotPri[i*pri[j]] = true;
            pw[i*pri[j]] = pw[i]*pw[pri[j]]%P;
            if (!(i%pri[j])) break;
        }
    }
    for (int i = 1; i <= n; i++) fac[i] = fac[i-1]*i%P;
    for (int i = 2; i <= n; i++) inv[i] = (P-P/i*inv[P%i]%P)%P;
    for (int i = 2; i <= n; i++) inv[i] = inv[i-1]*inv[i]%P;
}
lnt C(int n, int m) {return fac[n]*inv[m]%P*inv[n-m]%P;}
void get_Pnt_Val() {
    lnt k[MAX_M+5], b[MAX_M+5], invm;
    k[0] = 1, b[0] = 0, invm = Pow(m, P-2);
    for (int i = 1; i <= m+1; i++)
        k[i] = k[i-1]*invm%P, 
        b[i] = (b[i-1]*invm%P+pw[i])%P;
    int f = (m&1) ? -1 : 1; k[0] *= f, f = -f;
    for (int i = 1; i <= m+1; i++, f = -f)
        k[0] = (k[0]+f*C(m+1, i)*k[i]%P)%P, 
        b[0] = (b[0]+f*C(m+1, i)*b[i]%P)%P;
    F[0] = -b[0]*Pow(k[0], P-2)%P;
    for (int i = 1; i <= m; i++)
        F[i] = (k[i]*F[0]%P+b[i])%P;
}
#define getL(i) (i < m ? pre[m-i-1] : 1)
#define getR(i) (i > 0 ? suc[m-i+1] : 1)
lnt Poly_Inter() {
    lnt ret = 0;
    lnt pre[MAX_M+5], suc[MAX_M+5]; pre[0] = n-m, suc[m] = n;
    for (int i = 1; i <= m; i++) pre[i] = pre[i-1]*(n-m+i)%P;
    for (int i = m-1; i >= 1; i--) suc[i] = suc[i+1]*(n-m+i)%P;
    for (int i = 0, f = (m&1) ? -1 : 1; i <= m; i++, f = -f)
        ret = (ret+f*getL(i)*getR(i)%P*inv[i]%P*inv[m-i]%P*F[i]%P)%P;
    return ((Pow(m, n)*ret%P-F[0])%P+P)%P;
}
int main() {
    read(n), read(m), init(m+1);
    if (m == 1)
        printf("%lld\n", 1LL*n*(n+1)/2%P);
    else if (n <= m) {
        lnt pwm = m, sum = 0;
        for (int i = 1; i <= n; i++, pwm = pwm*m%P)
            sum = (sum+pw[i]*pwm%P)%P;
        printf("%lld\n", sum);
    } else
        get_Pnt_Val(), printf("%lld\n", Poly_Inter());
    return 0;
}