自然数幂和题解报告
原文链接：https://www.cnblogs.com/zsxuan/p/18016759
1. 线性逆元板子：https://www.luogu.com.cn/problem/P3811
题意：
线性求出 \(1 \sim n\) 模 \(m\) 的逆元。
\(1 \leq n \leq 3 \cdot 10^6, n < p < 2 \cdot 10^7\) 。
题解：
\(p\) 保证为质数且 \(p > n\) 意味着 \(gcd(p, n) = 1\) ，于是可以由裴蜀定理证明有逆元。
线性逆元
\(i\) 可由 \(p \bmod i\) 递推，当 \(i \geq 1\) 。于是 \(\bmod p\) 意义下：
且 \(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 \leq k \leq 2000, 1 \leq k,p \leq 10^9\) 。
题解 1 ：
算两次的差分思想：
\(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)\) 。
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\) 是完全积性函数，满足
于是借助线性筛有 \(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)\) 。
\(0 ~ n\) 的连续拉格朗日插值
求前缀和即得到 \(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\) 阶多项式
给出 \(n\) 个点，求 \(f(k) \bmod p\) 。
\(1 \leq 2 \cdot n \leq 10^3, 1 \leq x_i,y_i,k < 998244353\) 。
题解：
\(998244353\) 为质数，典型的存在逆元。
已知 \(n\) 阶多项式的拉格朗日插值
此时逆元无法线性预处理，需要使用快速幂。
将 \(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 \leq n \leq 1.6 \times 10^5 < m \leq 10^{8}\) 。
题解：
首先不难有连续插值公式
可以通过 \(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 \leq m \leq 50, 1 \leq a_i < n \leq 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}\) 点伤害，于是有一个差分的处理：
显然处理 \(a_0 = 0\) 。
于是击杀所有怪物可以获得的贡献为
后半部分 \(\sum_{i = 1}^{m + 1} \sum_{j = i}^{m} (a_j - a_{i - 1})^{m + 1}\) 可以直接 \(O(m^2)\) 计算，于是只需计算
不难发现 \(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 \leq 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\) 。
求 $ (g(a) + g(a + d) + g(a + 2d) + \cdots + g(a + nd)) \bmod p $ 。
\(1 \leq k \leq 100, 0 \leq a, n, d \leq 10^9\) ，\(p = 1234567891\) 。
题解：
经过验算 \(p = 1234567891\) 为一个 \(> 10^9\) 的质数。典型地存在逆元。
并且注意到 \(p + p > 2^{31} - 1\) 。
不妨试着直接写出式子
观察到可以写成一个式子
前置定理
幂和阶定理： \(\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\) 开始将贡献 \(1 n^0\) ，低阶贡献不影响阶数。
\(\square\)
推论： \(\sum_{i = 0}^{n} (a + bi)^k\) 是关于 \(n\) 的 \(k + 1\) 阶多项式 \(f(n) =\sum_{i = 0}^{k + 1} a_i n^{i}\) 。
证明：
已知多项式的阶数只有最高项的幂次确定，于是多项式乘以常数阶数不变。
由幂和阶定理，原式结果为关于 \(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\) 阶多项式。
\(\square\)
多项式等差前缀和升阶定理：关于 \(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\) 为任意非零整数。
证明：
由幂和阶定理的推论，原式为关于 \(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\) ，低阶贡献不影响阶数。
\(\square\)
由这几个典型定理，还能观察到
\(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)\) 。
连续拉格朗日插值
其中 \(\frac{\prod_{j = 0}^{n} (x - j)}{x - i}\) 求前后缀积即可 \(O(1)\) 预处理。
特别注意这里的 \(x\) 可能达到 \(long\ long\) ，可以先乘以 \(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 \leq N, R_i, M \leq 100, 1 \leq U_i \leq 10^{9}\)
题解：
数据范围似乎可以接受 DP ，于是尝试考虑 DP 做法。
设 \(f_{i, j}\) 为考虑了 \(i\) 门必修课，有 \(j\) 人被 B 神碾压的方案数。
考虑 \(f_{i, j}\) 如何转移向 \(i + 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 神。 -
考虑第 \(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\) 种可能的分数。
\(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})\) 。
于是得到转移方程
注意整理：
- 一个显然的枚举上界至少是 \(j + k \leq N - 1\) ，即（\(j + k \leq N - 1 \Rightarrow k \leq N - 1 - j\)）
- 原本被碾压的 \(k\) 个人现在不再被碾压，则这 \(k\) 个人至少需要满足的约束是这门课的成绩高于 B 神，即 \(k \leq R_{i} - 1\) 。
- 现在能被碾压的人数严格上界是 \(N - R_{i}\) ，即 \(j \leq 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(N M)\) 的复杂度预处理出 \(\forall 1 \leq i \leq 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;
}
