自然数幂和题解报告

原文链接：https://www.cnblogs.com/zsxuan/p/18016759

1. 线性逆元板子：https://www.luogu.com.cn/problem/P3811

题意：
线性求出 $1\sim n$ 模 $m$ 的逆元。
$1\le n\le 3\cdot {10}^6,n<p<2\cdot {10}^7$ 。

题解：
$p$ 保证为质数且 $p>n$ 意味着 $gcd(p,n)=1$ ，于是可以由裴蜀定理证明有逆元。

线性逆元
$i$ 可由 $pmodi$ 递推，当 $i\ge 1$ 。于是 $modp$ 意义下：

$$\begin{array}{rl}\lfloor \frac{p}{i}\rfloor i+pmodi& \equiv 0\\ i^{-1}& \equiv (p-\lfloor \frac{p}{i}\rfloor p)(pmodi{)}^{-1}\end{array}$$

且 $1^{-1}=1$ 。

时间复杂度 $O(n)$ 。

view

#include <bits/stdc++.h>
 
typedef long long ll;
 
const int N = 3E6 + 5;
int MOD;
 
int add(int a, int b) { return ( 1LL * (1LL * a + b) % MOD + MOD) % MOD; }
int mul(int a, int b) { return (1LL * a * b % MOD + MOD) % MOD; }
 
int inv[N], fac[N], ifac[N];
struct invInit {
    invInit() {
        fac[0] = ifac[0] = inv[1] = 1;
        for (int i = 1; i < N; i++) {
            if (i > 1) inv[i] = 1LL * (MOD - MOD / i) * inv[MOD % i] % MOD;
            fac[i] = 1LL * fac[i - 1] * i % MOD;
            ifac[i] = 1LL * ifac[i - 1] * inv[i] % MOD;
        }
    }
};
 
signed main() {
#ifndef ONLINE_JUDGE
freopen("IO/in", "r", stdin); freopen("IO/out", "w", stdout);
#endif 
    int n, p; std::cin >> n >> p;
    MOD = p;
    invInit __init__inv;
    for (int i = 1; i <= n; i++) {
        std::cout << inv[i] << "\n";
    }
    return 0;
}
 

2. 自然数幂和 1 ：http://oj.daimayuan.top/course/22/problem/1106

题意：
求 $\sum _{i=1}^n{i}^k$ 。
$1\le k\le 2000,1\le k,p\le {10}^9$ 。

题解 1 ：
算两次的差分思想：

$$\begin{array}{rl}(n+1{)}^{k+1}-n^{k+1}& =\sum _{i=0}^{k+1}(\genfrac{}{}{0ex}{}{k+1}{i})n^i-n^{k+1}=\sum _{i=0}^k(\genfrac{}{}{0ex}{}{k+1}{i})n^i\\ n^{k+1}-(n-1{)}^{k+1}& =\sum _{i=0}^k(\genfrac{}{}{0ex}{}{k+1}{i})(n-1{)}^i\\ (n-1{)}^{k+1}-(n-2{)}^{k+1}& =\sum _{i=0}^k(\genfrac{}{}{0ex}{}{k+1}{i})(n-2{)}^i\\ \cdots \\ 2^{k+1}-1^{k+1}& =\sum _{i=0}^k(\genfrac{}{}{0ex}{}{k+1}{i})1^i\\ \mathbf{\text{两边求和}}\\ (n+1{)}^{k+1}-1& =\sum _{j=1}^n\sum _{i=0}^k(\genfrac{}{}{0ex}{}{k+1}{i})j^i=\sum _{i=0}^k(\genfrac{}{}{0ex}{}{k+1}{i})S_i(n)\\ \mathbf{\text{分离出最高项}}\\ (n+1{)}^{k+1}-1& =\sum _{i=0}^{k-1}(\genfrac{}{}{0ex}{}{k+1}{i})S_i(n)+(k+1)S_k(n)\\ S_k(n)=\frac{(n+1{)}^{k+1}-1-\sum _{i=0}^{k-1}(\genfrac{}{}{0ex}{}{k+1}{i})S_i(n)}{k+1}\end{array}$$

$S_0(n)=n$ 。

时间复杂度显然的 $O(k^2)$ 。

view

#include <bits/stdc++.h>
 
typedef long long ll;
 
const int N = 2E3 + 10;
 
int k, n, p;
int s[N], pw[N];
 
int inv[N], fac[N], ifac[N];
struct invInit {
    invInit() {
        inv[1] = fac[0] = ifac[0] = 1;
        for (int i = 1; i < N; i++) {
            if (i > 1) inv[i] = 1LL * (p - p / i) * inv[p % i] % p;
            fac[i] = 1LL * fac[i - 1] * i % p;
            ifac[i] = 1LL * ifac[i - 1] * inv[i] % p;
        }
    }
};
 
int comb(int n, int m) {
    if (n < 0 || m < 0 || n < m) return 0;
    return 1LL * fac[n] * ifac[n - m] % p * ifac[m] % p;
}
 
signed main() {
#ifndef ONLINE_JUDGE
freopen("IO/in", "r", stdin); freopen("IO/out", "w", stdout);
#endif 
    std::cin >> k >> n >> p;
    invInit __init__inv;
    pw[0] = 1;
    for (int i = 1; i <= k + 1; i++) pw[i] = 1LL * pw[i - 1] * (n + 1) % p;
    s[0] = n;
    for (int i = 1; i <= k; i++) {
        int P = (pw[i + 1] - 1 + p) % p;
        for (int j = 0; j < i; j++)
            P = (1LL * P - 1LL * comb(i + 1, j) * s[j] % p + p) % p;
        int Q = inv[i + 1];
        s[i] = 1LL * P * Q % p;
    }
    std::cout << s[k] << "\n";
	return 0;
}

题解 2 ：
$\sum _{i=1}^n{i}^k$ 是关于 $n$ 的 $k+1$ 次多项式，由题解 1 的递推式可归纳。
于是至少可以 $O(k^2)$ 进行朴素拉格朗日插值。
不妨由递推式 $S_k(n+1)=S_k(n)+(n+1{)}^k$ 直接递推求出 $k+2$ 个点对。时间复杂度 $O(k\log k)$ ，其中 $\log k$ 由快速幂贡献。
复杂度瓶颈在 $O(k^2)$ ，再贡献一个 $O(\log p)$ 的逆元。于是总复杂度为 $O(k^2\log p)$ 。

$$\begin{array}{r}f(x)=\sum _{i=0}^n{y}_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}\end{array}$$

view

#include <bits/stdc++.h>
 
typedef long long ll;
 
const int N = 2010;
int MOD;
 
int add(int a, int b) { return ( 1LL * (1LL * a + b) % MOD + MOD) % MOD; }
int mul(int a, int b) { return (1LL * a * b % MOD + MOD) % MOD; }
int qpow(int a, ll n) { int res=1; for(;n;a=mul(a,a),n>>=1) if(n&1) res=mul(res,a); return res; }
 
int k, n, p;
int x[N], y[N];
 
int lagrange(int k, int * x, int *y, int n) {
    int fk = 0;
    for (int i = 0; i <= n; i ++) {
        int P = y[i], Q = 1;
        for (int j = 0; j <= n; j++) if (j != i) {
            P = mul(P, add(k, -x[j]));
            Q = mul(Q, add(x[i], -x[j]));
        }
        Q = qpow(Q, MOD - 2);
        fk = add(fk, mul(P, Q));
    }
    return fk;
}
 
signed main() {
#ifndef ONLINE_JUDGE
freopen("IO/in", "r", stdin); freopen("IO/out", "w", stdout);
#endif 
    std::cin >> k >> n >> p;
    MOD = p;
    for (int i = 0; i <= k + 1; i++) x[i] = i;
    // S_{n + 1}(k) = S_{n}(k) + (n + 1)^k, S_{0}(k) = 0
    y[0] = 0;
    for (int i = 1; i <= k + 1; i++) y[i] = add(y[i - 1], qpow(i, k));
    int ans = lagrange(n, x, y, k + 1);
    std::cout << ans << "\n";
    
	return 0;
}

3. 自然数幂和 2 ：https://www.luogu.com.cn/problem/CF622F

续借上文的自然数幂和 1 ，考虑如何优化。

首先是 $k+2$ 个点对，直接递推需要使用快速幂，存在其他方法吗？
观察到 $f(i)=i^k$ 是完全积性函数，满足

$$\begin{array}{r}\{\begin{array}{ll}f(1)& =1^k\\ f(ab)& =f(ab)\end{array}\end{array}$$

于是借助线性筛有 $f(p_1\times i)=f(p_1)\times f(i)$ ，在合数时 $O(1)$ 递推，质数时 $O(\log n)$ 快速幂计算。
显然筛法的瓶颈在于质数密度。根据质数定理 $n$ 以内的质数个数约为 $\frac{n}{\ln n}$ ，于是筛法时间复杂度为 $T({\log }_{2}n\frac{n}{\ln n}=n\frac{{\log }_{2}n}{\ln n}=n\frac{\ln n/\ln 2}{\ln n}=n/\ln 2)=O(n)$ 。

$0n$ 的连续拉格朗日插值

$$f(x)=\sum _{i=0}^n\frac{y_i(-1{)}^{n-i}}{i!(n-i)!}\frac{\prod _{j=0}^n(x-j)}{x-i}$$

求前缀和即得到 $k+2$ 个点对。

view

#include <bits/stdc++.h>
 
typedef long long ll;
 
const int N = 1E6 + 10;
const int MOD = 1E9 + 7;
 
int n, k;
int y[N], inv[N], fac[N], ifac[N];
int p[N], pr[N], tot;
 
int add(int a, int b) { return ( 1LL * (1LL * a + b) % MOD + MOD) % MOD; }
int mul(int a, int b) { return (1LL * a * b % MOD + MOD) % MOD; }
int qpow(int a, ll n) { int res=1; for(;n;a=mul(a,a),n>>=1) if(n&1) res=mul(res,a); return res; }
 
struct eulaSieveInit {
    eulaSieveInit() {
        y[1] = 1;
        for (int i = 2; i < N; i++) {
            if (!p[i]) p[i] = i, pr[++tot] = i, y[i] = qpow(i, k);
            for (int j = 1; j <= tot && i * pr[j] < N; j++) {
                p[i * pr[j]] = pr[j];
                y[i * pr[j]] = mul(y[i], y[pr[j]]);
                if (p[i] == pr[j]) {
                    break;
                }
            }
        }
    }
};
 
struct invInit {
    invInit () {
        inv[1] = fac[0] = ifac[0] = 1;
        for (int i = 1; i < N; i++) {
            if (i > 1) inv[i] = mul(MOD - MOD / i, inv[MOD % i]);
            fac[i] = mul(fac[i - 1], i);
            ifac[i] = mul(ifac[i - 1], inv[i]);
        }
    }
};
 
int consecutive_lagrange(int k, int y[], int n) {
    int fk = 0;
    std::vector<int> pre(n + 2), suf(n + 2);
    pre[0] = k; suf[n + 1] = 1;
    for (int i = 1; i <= n; i++) pre[i] = mul( pre[i - 1] , add( k , -i ) );
    for (int i = n; i >= 0; --i) suf[i] = mul( suf[i + 1] , add( k , -i ) );
    for (int i = 0; i <= n; i++) {
        int sgn = ((n - i) & 1) ? -1 : 1;
        int res = mul( y[i] , sgn );
        int P = mul( i > 0 ? pre[i - 1] : 1, suf[i + 1] );
        int Q = mul(ifac[i] , ifac[n - i]);
        res = mul( res , mul( P , Q ) );
        fk = add( fk , res );
    }
    return fk;
}
 
signed main() {
#ifndef ONLINE_JUDGE
freopen("IO/in", "r", stdin); freopen("IO/out", "w", stdout);
#endif 
    std::cin >> n >> k;
    invInit __init__inv;
    y[0] = 0;
    eulaSieveInit __init_eulaSieve;
    for (int i = 1; i <= k + 1; i++) y[i] = add(y[i], y[i - 1]);
    int ans = consecutive_lagrange(n, y, k + 1);
    std::cout << ans << "\n";
    return 0;
}

4. 拉格朗日插值板子 1 ：https://www.luogu.com.cn/problem/P4781

题意：
对于一个关于 $x$ 的 $n$ 阶多项式

$$f(x)=a_0{x}_0+a_1{x}_1{a}_2{x}_2+\cdots +a_n{x}_n$$

给出 $n$ 个点，求 $f(k)modp$ 。

$1\le 2\cdot n\le {10}^3,1\le x_i,y_i,k<998244353$ 。

题解：

$998244353$ 为质数，典型的存在逆元。

已知 $n$ 阶多项式的拉格朗日插值

$$f(x)=\sum _{i=0}^n{y}_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}$$

此时逆元无法线性预处理，需要使用快速幂。

将 $k$ 带入 $n-1$ 阶多项式的插值公式即可。

时间复杂度 $O(n^2\log p)$ 。

view

#include <bits/stdc++.h>
 
typedef long long ll;
 
const int N = 2010;
const int MOD = 998244353;
 
int add(int a, int b) { return ( 1LL * (1LL * a + b) % MOD + MOD) % MOD; }
int mul(int a, int b) { return (1LL * a * b % MOD + MOD) % MOD; }
int qpow(int a, ll n) { int res=1; for(;n;a=mul(a,a),n>>=1) if(n&1) res=mul(res,a); return res; }
 
int k, n;
int x[N], y[N];
 
int lagrange(int k, int * x, int *y, int n) {
    int fk = 0;
    for (int i = 0; i <= n; i ++) {
        int P = y[i], Q = 1;
        for (int j = 0; j <= n; j++) if (j != i) {
            P = mul(P, add(k, -x[j]));
            Q = mul(Q, add(x[i], -x[j]));
        }
        Q = qpow(Q, MOD - 2);
        fk = add(fk, mul(P, Q));
    }
    return fk;
}
 
signed main() {
#ifndef ONLINE_JUDGE
freopen("IO/in", "r", stdin); freopen("IO/out", "w", stdout);
#endif 
    std::cin >> n >> k;
    for (int i = 0; i < n; i++) {
        std::cin >> x[i] >> y[i];
    }
    std::cout << lagrange(k, x, y, n - 1) << "\n";
    return 0;
}
 

5. 拉格朗日插值板子 2 ：https://www.luogu.com.cn/problem/P5667

题意：
给定一个 $n$ 阶多项式的 $n+1$ 个点值 $(0,f(0),(1,f(1),(2,f(2),\cdots ,(n,f(n))$ 。求 $n+1$ 个值 $f(m),f(m+1),f(m+2),\cdots ,f(m+n)$ 。答案对 $998244353$ 取模。
$1\le n\le 1.6\times {10}^5<m\le {10}^8$ 。

题解：
首先不难有连续插值公式

$$f(x)=\sum _{i=0}^n{y}_i(-1{)}^{n-i}\frac{\prod _{j=0}^n(x-j)}{x-i}\frac{1}{i!(n-i)!}$$

可以通过 $O(n)$ 复杂度插值出 $f(m)$ 。

但求出 $f(m)\sim f(m+n)$ 将会是 $O(n^2)$ 的时间复杂度，不可接受。

6. 拉格朗日插值水题 1 ：https://www.luogu.com.cn/problem/P4593

题意：
有一些怪物，每只怪物血量各不相同。血量最高的怪物为 $n$ ，有 $m$ 个血量 $a_1,a_2,a_2,\cdots ,a_m$ 代表这些血量没有出现过。即此外 $\{1,2,3,\cdots ,n\}{a}_i$ 的血量的怪物各有一只。怪物血量为 $0$ 代表死亡。

现在手中有无限张“亵渎”卡，一张卡的效果为
1.对所有怪物造成一点伤害
2.每有一只怪物死亡，重复一次法术

每次使用“亵渎”卡后，每只收到伤害怪物会贡献的分数为 $x^k$ 。$x$ 使用“亵渎”前这只怪物的血量，$k$ 为击杀所有怪物需要的“亵渎”张数。

$1\le m\le 50,1\le a_i<n\le {10}^{13}$ 。

题解：
首先确定 $k$ 为常数，考虑如何计算 $k$ 。

让 $a$ 数组为升序。

已知怪物血量不重复。观察空缺血量 $a_i$ ，若 $a_i-1$ 为空缺血量，则 $a_i$ 需要消耗一张“亵渎。”若 $a_i-1$ 为存在血量，则 $a_i$ 至 $a_i-1$ 起的一段存在血量的前缀需要消耗一张“亵渎”。

于是直到最后一个空缺血量 $a_m$ ，共需消耗 $m$ 张“亵渎”。考虑最后一段存在血量的后缀，可以得出 $m+1$ 张“亵渎”可以击杀所有怪物，即 $k=m+1$ 。

任意一只怪物对每次“亵渎”卡的贡献的 $x$ 不同。任意一次“亵渎”卡使用时每只怪物的 $x$ 可能不同。

不难发现死亡消失的怪物不会产生贡献。

单独考虑第 $i$ 张"亵渎"使用前的情况，不难观察到所有存货的怪物都受过 $a_{i-1}$ 点伤害，于是有一个差分的处理：

$$\begin{array}{rl}& \sum _{j=1}^n(j-a_{i-1}{)}^{m+1}-\sum _{j=1}^m(a_j-a_{i-1}{)}^{m+1}\\ \mathbf{\text{舍去非正血量下标的贡献}}\\ & =\sum _{j=a_i+1}^n(j-a_{i-1}{)}^{m+1}-\sum _{j=i}^m(a_j-a_{i-1}{)}^{m+1}\\ & =\sum _{j=1}^{n-a_{i-1}}{j}^{m+1}-\sum _{j=i}^m(a_j-a_{i-1}{)}^{m+1}\end{array}$$

显然处理 $a_0=0$ 。

于是击杀所有怪物可以获得的贡献为

$$\sum _{i=1}^{m+1}(\sum _{j=1}^{n-a_{i-1}}{j}^{m+1}-\sum _{j=i}^m(a_j-a_{i-1}{)}^{m+1})$$

后半部分 $\sum _{i=1}^{m+1}\sum _{j=i}^m(a_j-a_{i-1}{)}^{m+1}$ 可以直接 $O(m^2)$ 计算，于是只需计算

$$\sum _{i=1}^{m+1}\sum _{j=1}^{n-a_{i-1}}{j}^{m+1}$$

不难发现 $f(n-a_{i-1})=\sum _{j=1}^{n-a_{i-1}}{j}^{m+1}$ 是 $m+2$ 阶多项式，$g(m+1)=\sum _{i=1}^{m+1}\sum _{j=1}^{n-a_{i-1}}{j}^{m+1}$ 是 $m+3$ 阶多项式。

瓶颈在于 $n\le {10}^{13}$ ，于是可以线性连续插值 $O(m)$ 得到 $f$ ，这个处理是经典的线性筛+连续拉格朗日插值。再直接前缀和以 $O(m^2)$ 算出 $g(m+1)$ 。

注意线性筛的写法，每次需要初始化质数数组，不然只能筛一次函数。（被演了）

总复杂度 $O(m+1)$ 。

view

#include <bits/stdc++.h>
 
typedef long long ll;
 
const int N = 60;
const int MOD = 1E9 + 7;
 
ll n, m;
int a[N], y[N];
int p[N], pr[N], tot;
int inv[N], fac[N], ifac[N];
 
int add(int a, int b) { return ( 1LL * (1LL * a + b) % MOD + MOD) % MOD; }
int mul(int a, int b) { return (1LL * a * b % MOD + MOD) % MOD; }
int qpow(int a, ll n) { int res=1; for(;n;a=mul(a,a),n>>=1) if(n&1) res=mul(res,a); return res; }
 
struct eulaSieveInit {
    eulaSieveInit() {
        tot = 0;
        for (int i = 1; i < N; i++) p[i] = 0;
        y[1] = 1;
        for (int i = 2; i < N; i++) {
            if (!p[i]) p[i] = i, pr[++tot] = i, y[i] = qpow(i, m + 1);
            for (int j = 1; j <= tot && i * pr[j] < N; j++) {
                p[i * pr[j]] = pr[j];
                y[i * pr[j]] = mul(y[i], y[pr[j]]);
                if (p[i] == pr[j]) {
                    break;
                }
            }
        }
    }
};
 
struct invInit {
    invInit() {
        inv[1] = fac[0] = ifac[0] = 1;
        for (int i = 1; i < N; i++) {
            if (i > 1) inv[i] = mul(MOD - MOD / i, inv[MOD % i]);
            fac[i] = mul(fac[i - 1], i);
            ifac[i] = mul(ifac[i - 1], inv[i]);
        }
    }
};
 
int consecutive_lagrange(int k, int *y, int n) {
    int fk = 0;
    std::vector<int> pre(n + 2), suf(n + 2);
    for (int i = 0; i <= n; i++) pre[i] = i > 0 ? mul(pre[i - 1], k - i) : k;
    suf[n + 1] = 1; for (int i = n; i >= 0; --i) suf[i] = mul(suf[i + 1], k - i);
    for (int i = 0; i <= n; i++) {
        int sgn = ((n - i) & 1) ? -1 : 1;
        int res = mul( i > 0 ? pre[i - 1] : 1, suf[i + 1]);
        int P = mul(y[i], sgn);
        int Q = mul(ifac[i], ifac[n - i]);
        fk = add(fk, mul(res, mul(P, Q)));
    }
    return fk;
}
 
signed main() {
#ifndef ONLINE_JUDGE
freopen("IO/in", "r", stdin); freopen("IO/out", "w", stdout);
#endif 
    int tc = 0; std::cin >> tc;
    while (tc--) {
        std::cin >> n >> m;
        invInit __init__inv;
        eulaSieveInit __init__eulaSieve;
        for (int i = 1; i < N; i++) y[i] = add(y[i], y[i - 1]);
        for (int i = 1; i <= m; i++) {
            std::cin >> a[i];
        }
        std::sort(a + 1, a + m + 1);
        int ans = 0;
        for (int i = 1; i <= m + 1; i++) {
            ans = add(ans, consecutive_lagrange(add(n , - a[i - 1]), y, m + 2));
        }
        for (int i = 1; i <= m + 1; i++) {
            for (int j = i; j <= m; j++) {
                ans = add(ans, -qpow(add(a[j], - a[i - 1]), m + 1));
            }
        }
        std::cout << ans << "\n";
    }
 
    return 0;
}

7. 拉格朗日插值水题 2 ：https://vjudge.net/problem/BZOJ-3453

题意：
给定 $k,a,n,d,p$ 。

$$\begin{array}{rl}f(i)& =1^k+2^k+3^k+\cdots +i^k\\ g(x)& =f(1)+f(2)+f(3)+\cdots +f(x)\end{array}$$

求 $(g(a)+g(a+d)+g(a+2d)+\cdots +g(a+nd))modp$ 。

$1\le k\le 100,0\le a,n,d\le {10}^9$ ，$p=1234567891$ 。

题解：

经过验算 $p=1234567891$ 为一个 $>{10}^9$ 的质数。典型地存在逆元。

并且注意到 $p+p>2^{31}-1$ 。

不妨试着直接写出式子

$$\begin{array}{rl}h(n)& =\sum _{i=0}^{n}g(a+id)\\ g(m)& =\sum _{i=1}^{m}f(i)\\ f(l)& =\sum _{i=1}^l{i}^k\end{array}$$

观察到可以写成一个式子

$$\begin{array}{rl}h(n)& =\sum _{i=0}^n\sum _{j=1}^{a+id}f(j)\\ & =\sum _{i=0}^n\sum _{j=1}^{a+id}\sum _{l=1}^j{l}^k\end{array}$$

前置定理

幂和阶定理： $\sum _{i=b}^n{i}^{k}s.t.b=0,1$ 是关于 $n$ 的 $k+1$ 阶多项式。
证明：
由积累（大雾），通过 $\sum _{i=1}^n{i}^k$ 的组合递推式可归纳出这是一个关于 $n$ 的 $k+1$ 次多项式 $f(n)=\sum _{i=0}^{k+1}{a}_i{n}^i$ 。
递推式计算过程略。
若 $i$ 从 $0$ 开始将贡献 $1n^0$ ，低阶贡献不影响阶数。
$◻$
推论： $\sum _{i=0}^n(a+bi{)}^k$ 是关于 $n$ 的 $k+1$ 阶多项式 $f(n)=\sum _{i=0}^{k+1}{a}_i{n}^i$ 。
证明：

$$\begin{array}{rl}\sum _{i=1}^n(u+vi{)}^k& =\sum _{i=1}^n\sum _{j=0}^k(\genfrac{}{}{0ex}{}{k}{j})u^{k-j}{v}^j{i}^j\\ & =\sum _{j=0}^k(\genfrac{}{}{0ex}{}{n}{k-j})u^{k-j}{v}^j\sum _{i=1}^n{i}^j\end{array}$$

已知多项式的阶数只有最高项的幂次确定，于是多项式乘以常数阶数不变。
由幂和阶定理，原式结果为关于 $n$ 的 $j+1$ 阶多项式 $\sum _{i=0}^{j+1}{a}_i{n}^i$ 的前缀和。只看最高阶多项式，为关于 $n$ 的 $k+1$ 阶多项式 $\sum _{i=0}^{k+1}{a}_i{n}^i$ 。
于是 $\sum _{i=1}^n(u+vi{)}^k$ 是关于 $n$ 的 $k+1$ 阶多项式。
$◻$
多项式等差前缀和升阶定理：关于 $u+vk$ 的 $n$ 阶多项式的 $m$ 阶等差前缀和 $\sum _{i=b}^{m}f(u+vi)s.t.b=0,1$ 是关于 $m$ 的 $n+1$ 阶多项式 $\sum _{i=0}^{n+1}{a}_i{m}^i$ 。$u$ 为任意整数，$v$ 为任意非零整数。
证明：

$$\begin{array}{r}\sum _{k=1}^m\sum _{i=0}^n{a}_i(u+vk{)}^i=\sum _{i=0}^n{a}_i\sum _{k=1}^m(u+vk{)}^i\end{array}$$

由幂和阶定理的推论，原式为关于 $m$ 的 $i+1$ 次多项式 $\sum _{j=0}^{i+1}{a}_j{m}^j$ 的前缀和。只看最高阶多项式，为关于 $m$ 的 $n+1$ 阶多项式 $\sum _{j=0}^{n+1}{a}_j{m}^j$ 。
于是 $\sum _{k=1}^m\sum _{i=0}^n{a}_i(u+vk{)}^i$ 是关于 $m$ 的 $n+1$ 阶多项式。
前缀和从 $0$ 开始将贡献 $C{m}^0$ ，低阶贡献不影响阶数。
$◻$

由这几个典型定理，还能观察到

$h(n)=\sum _{i=0}^n\sum _{j=1}^{a+id}\sum _{l=1}^j{l}^k$ 是关于 $n$ 的 $k+3$ 阶多项式。

我们可以考虑求出 $(0,h(0)),(1,h(1)),(2,h(2)),\cdots ,(k+4,h(k+4))$ 直接求出 $k+4$ 个 $h$ 的点对。

似乎有点难，但我们可以先算出差分数组。

耿直地，可以以低于 $O(k\log p)$ 求出 $k+2$ 个 $f$ 的连续点对，并 $O(k)$ 插值得到 $f(x)$ 。
通过 $O(k^2)$ 求出 $k+3$ 个 $g$ 的连续点对，并 $O(k)$ 插值得到 $g(x)$ 。
于是可通过 $O(k^2)$ 求出 $k+4$ 个 $h$ 的连续点对。

瓶颈是 $O(k^2)$ 的，$f$ 不需要线性筛点对。并且，不妨全部求 $k+4$ 个点对。

考虑简化，冷静分析。$h$ 是 $g$ 的等差前缀和，$g$ 是 $f$ 的前缀和。且前缀和不是瓶颈，瓶颈是等差前缀和。

考虑 $g(n)=\sum _{i=1}^{n}f(i)$ 。可以通过 $O(k\log p)$ 递推求出 $k+4$ 个 $f$ 的连续点对，然后 $O(k)$ 简单前缀和即可得到 $k+4$ 个 $g$ 的连续点对。

于是可以直接 $O(k)$ 插值出 $n+2$ 阶多项式 $g$ 。

继而有 $h(n)=\sum _{i=0}^{n}g(a+id)$ ，通过 $g$ 的等差前缀和求出 $O(k)$ 个 $h$ 点对，时间复杂度 $O(k^2)$ 。

再次进行 $O(k)$ 插值求出 $n+3$ 阶多项式 $h(x)$ 。

连续拉格朗日插值

$$\begin{array}{rl}f(x)& =\sum _{i=0}^n{y}_i\prod _{j\ne i}\frac{x-j}{i-j}\\ & =\sum _{i=0}^n{y}_i\frac{\prod _{j=0}^n(x-j)}{\prod _{j\ne i}(i-j)(x-i)}\\ & =\sum _{i=0}^n\frac{y_i(-1{)}^{n-i}}{i!\times (n-i)!}\frac{\prod _{j=0}^n(x-j)}{(x-i)}\end{array}$$

其中 $\frac{\prod _{j=0}^n(x-j)}{x-i}$ 求前后缀积即可 $O(1)$ 预处理。

特别注意这里的 $x$ 可能达到 $longlong$ ，可以先乘以 $1$ 并取模。
要特别注意 $x$ 的范围。

view

#include <bits/stdc++.h>
typedef long long ll;
 
const int N = 2E2 + 10;
const int MOD = 1234567891;
 
int k, a, n, d;
int inv[N], fac[N], ifac[N];
int f[N], g[N], h[N];
int add(int a, int b) { return ( 1LL * (1LL * a + b) % MOD + MOD) % MOD; }
int mul(int a, int b) { return (1LL * a * b % MOD + MOD) % MOD; }
int qpow(int a, ll n) { int res = 1; for (; n; a=mul(a, a), n>>=1) if(n&1) res=mul(res, a); return res; }
 
struct invInit {
    invInit() {
        fac[0] = ifac[0] = inv[1] = 1;
        for (int i = 1; i < N; i++) {
            if (i > 1) inv[i] = 1LL * (MOD - MOD / i) * inv[MOD % i] % MOD;
            fac[i] = 1LL * fac[i - 1] * i % MOD;
            ifac[i] = 1LL * ifac[i - 1] * inv[i] % MOD;
        }
    }
};
 
inline int consecutive_lagrange(int k, int *y, int n){
	std::vector<int> pre(n + 2), suf(n + 2);
    pre[0] = k; suf[n + 1] = 1;
    for (int i = 1; i <= n; i++) { pre[i] = mul(pre[i - 1] , add( k , -i )); }
    for (int i = n; i >= 0; --i) { suf[i] = mul(suf[i + 1], add( k, -i) ); }
    int fk = 0;
	for(int i = 0; i <= n; i++){
        int sgn = ((n - i) & 1) ? -1 : 1;
        int res = mul(i > 0 ? pre[i - 1] : 1, suf[i + 1]);
        int P = mul( sgn , y[i] );
        int Q = mul(ifac[i] , ifac[n - i]);
        fk = add( fk , mul( res, mul( P , Q ) ) );
	}
    return fk;
}
 
signed main() {
#ifndef ONLINE_JUDGE
freopen("IO/in", "r", stdin); freopen("IO/out", "w", stdout);
#endif
    invInit init;
    int tc = 0; std::cin >> tc;
    while (tc--) {
        std::cin >> k >> a >> n >> d;
        
        // // f(n + 1) = f(n) + (n + 1)^k
        f[0] = 0; // k > 0
        // // g(n) = sum_{i = 1}^{n} f(i)
        g[0] = 0;
        for (int i = 1; i <= k + 3; i++) f[i] = add( f[i - 1] , qpow(i, k) );
        for (int i = 1; i <= k + 3; i++) g[i] = add( g[i - 1] , f[i] );
        // // h(n) = sum_{i = 0}^{n} g(a + di)
        h[0] = consecutive_lagrange( mul(a, 1) , g, k + 2 );
        for (int i = 1; i <= k + 3; i++) {
            h[i] = add( h[i - 1], consecutive_lagrange( add( a, mul( i , d ) ) , g, k + 2 ) );
        }
        int h_n = consecutive_lagrange(n, h, k + 3);
        std::cout << h_n << "\n";
    }
    
    return 0;
}
 

8. 拉格朗日插值水题 3 ：https://www.luogu.com.cn/problem/P3270

题意：G 班有 $N$ 名同学，$M$ 门必修课。有 $K$ 位同学被 B 神碾压。一个同学若被 B 神碾压，则他的所有必修课分数均小于等于 B 神的必修课分数，但 B 神不碾压自己。现在给出 $M$ 门必修课的最高分 $U_i$ 和 B 神的课程排名 $R_i$ 。需要注意第 $i$ 门必修课中,严格有 $N-R_i$ 个人分数“大于” B 神，$R_i-1$ 个人分数“小于等于” B 神。
询问 G 班有多少种可能的分数情况。

$1\le N,R_i,M\le 100,1\le U_i\le {10}^9$

题解：
数据范围似乎可以接受 DP ，于是尝试考虑 DP 做法。
设 $f_{i,j}$ 为考虑了 $i$ 门必修课，有 $j$ 人被 B 神碾压的方案数。

考虑 $f_{i,j}$ 如何转移向 $i+1$ 。

1. 考虑排名的分布情况。
   若考虑了 $i+1$ 门必修课，原 $j$ 个人在前 $i$ 门必修课被 B 神碾压的情况下第 $i+1$ 门必修课可能超过 B 神，于是被 B 神碾压的人数可能变少。
   原来被 B 神碾压的 $j$ 个人有 $k$ 个人不再被碾压，于是从原来被碾压的 $j$ 个人中需要选 $k$ 个人，这门必修课成绩需要大于 B 神，另外的人这门必修课小于等于 B 神。
   此时第 $i+1$ 门必修课成绩高于 B 神的位置还剩 $R_{i+1}-1-k$ 个，于是从原来没被碾压的 $N-1-j$ 个人中选 $R_{i+1}-i-k$ 个人，这门必修课成绩需要大于 B 神，另外的人这门必修课小于 B 神。

2. 考虑第 $i+1$ 门课的排名确定后，每个排名上的人的分数分配情况。
   枚举 B 神第 $i+1$ 门课的可能成绩为 $j$ 。
   则在这门课中，有 $R_{i+1}-1$ 个排名高于 B 神的人分数范围在 $[j+1,U_{i+1}]$ ，每个人有 $U_{i+1}-j$ 种可能的分数。
   $N-R_{i+1}$ 个排名不高于 B 神的人分数范围在 [1, j] ，每个人有 $j$ 种可能的分数。

$$\begin{array}{r}{T}_{i+1}=\sum _{j=1}^{U_{i+1}}(U_{i+1}-j{)}^{R_{i+1}-1}{j}^{N-R_{i+1}}\end{array}$$

$U_{i+1}$ 的瓶颈很高，但仔细分析这一坨玩意是什么？
关于同一个变量的高阶多项式加上低阶多项式不影响阶数，不需要真的使用二项式定理展开。
比较直观地可以只考虑能取到的最高次 $j^{N-R_{i+1}}{j}^{R_{i+1}-1}=j^{N-1}$ 。这是一个经典的关于 $U_{i+1}$ 的 $N$ 阶多项式。
令 $T_{i+1}=f_N(U_{i+1})$ 。

于是得到转移方程

$$f_{i+1,j-k}=f_{i,j}(\genfrac{}{}{0ex}{}{j}{k})(\genfrac{}{}{0ex}{}{N-1-j}{R_{i+1}-1-k})T_{i+1}$$

注意整理：

1. 一个显然的枚举上界至少是 $j+k\le N-1$ ，即（$j+k\le N-1\Rightarrow k\le N-1-j$）
2. 原本被碾压的 $k$ 个人现在不再被碾压，则这 $k$ 个人至少需要满足的约束是这门课的成绩高于 B 神，即 $k\le R_i-1$ 。
3. 现在能被碾压的人数严格上界是 $N-R_i$ ，即 $j\le N-R_i$ 。

$$f_{i,j}=\sum _{k=0}^{min(N-1-j,N-R_i)}{f}_{i-1,j+k}(\genfrac{}{}{0ex}{}{j+k}{k})(\genfrac{}{}{0ex}{}{N-1-j-k}{R_i-1-k})T_{i}s.t.1\le i\le m,0\le j\le N-R_i$$

考虑 DP 数组的良序情况，第 $0$ 门必修课所有人都是 $0$ 分，B 神碾压所有人，即 $f_{0,n-1}=1$ 。

组合数可以 $O(N)$ 预处理后 $O(1)$ 计算。
$T_i=f_N(U_i+1)$ 可以 $O(N)$ 求出 $N+1$ 组点对，然后 $O(N)$ 插值出 $f_N(U_i+1)$ 。不同于 $DP$ 数组的良序处理，于是可以使用 $O(NM)$ 的复杂度预处理出 $\mathrm{\forall }1\le i\le M,T_i=f_N(U_i+1)$ ，然后 $O(1)$ 计算。
于是转移复杂度为 $O(1)$ ，状态复杂度为 $O(N^{2}M)$ 。故总时间复杂度为 $O(N^{2}M)$ 。

view

#include <bits/stdc++.h>
 
typedef long long ll;
 
const int N = 210;
const int CN = 110, CM = 110;
const int MOD = 1E9 + 7;
 
int add(int a, int b) { return ( 1LL * (1LL * a + b) % MOD + MOD) % MOD; }
int mul(int a, int b) { return (1LL * a * b % MOD + MOD) % MOD; }
int qpow(int a, ll n) { int res = 1; for (; n; a=mul(a, a), n>>=1) if(n&1) res=mul(res, a); return res; }
 
int U[CM], R[CM], f[CM][CN], T[CM];
int inv[N], fac[N], ifac[N], y[N];
int n, m, k;
 
struct invInit {
    invInit() {
        inv[1] = fac[0] = ifac[0] = 1;
        for (int i = 1; i < N; i++) {
            if (i > 1) inv[i] = mul(MOD - MOD / i, inv[MOD % i]);
            fac[i] = mul(fac[i - 1], i);
            ifac[i] = mul(ifac[i - 1], inv[i]);
        }
    }
};
 
int comb(int n, int m) {
    if (n < 0 || m < 0 || n < m) return 0;
    return mul(fac[n], mul(ifac[n - m], ifac[m]));
}
 
int consecutive_lagrange(int k, int *y, int n) {
    int fk = 0;
    std::vector<int> pre(n + 2), suf(n + 2);
    pre[0] = k; suf[n + 1] = 1;
    for (int i = 1; i <= n; i++) pre[i] = mul(pre[i - 1], k - i);
    for (int i = n; i >= 0; --i) suf[i] = mul(suf[i + 1], k - i);
    for (int i = 0; i <= n; i++) {
        int sgn = ((n - i) & 1) ? -1 : 1;
        int res = mul( i > 0 ? pre[i - 1] : 1, suf[i + 1]);
        int P = mul(y[i], sgn);
        int Q = mul(ifac[i], ifac[n - i]);
        fk = add(fk, mul(res, mul(P, Q)));
    }
    return fk;
}
 
signed main() {
#ifndef ONLINE_JUDGE
freopen("IO/in", "r", stdin); freopen("IO/out", "w", stdout);
#endif 
    invInit __init__inv;
    std::cin >> n >> m >> k;
    for (int i = 1; i <= m; i++) {
        std::cin >> U[i];
    }
    for (int i = 1; i <= m; i++) {
        std::cin >> R[i];
    }
 
    for (int i = 1; i <= m; i++) {
        y[0] = 0;
        for (int j = 1; j <= n; j++) {
            int res = mul(qpow(add(U[i], - j), add(R[i], - 1)), qpow(j, add(n , - R[i])));
            y[j] = add(y[j - 1], res); 
        }
        T[i] = consecutive_lagrange(U[i], y, n);
    }
    f[0][n - 1] = 1;
    for (int i = 1; i <= m; i++) {
        for (int j = 0; j <= n - R[i]; j++) {
            for (int l = 0; l <= std::min(n - 1 - j, R[i] - 1); l++) {
                int res = f[i - 1][j + l];
                res = mul(res, mul(comb(j + l, l), comb(n - 1 - j - l, R[i] - 1 - l)));
                res = mul(res, T[i]);
                f[i][j] = add(f[i][j], res);
            }
        }
    }
    std::cout << f[m][k] << "\n";
    return 0;
}