闲话 23.1.8
闲话
怎么了的 洛谷不能交太长的代码是吧
经过测试大概码长分界线在 11369 ~ 11381
而且发现 ABC284 的题已经整一天了都没上
洛谷是不是在作什么大动作？
今日推歌是《ぼくらの16bit戦争(Remaster ver)》sasakure.UK feat.GUMI
推荐博客：多项式 I：拉格朗日插值与快速傅里叶变换, Alex_Wei
拉格朗日插值例题和做法不是拉格朗日插值的题
感觉拉格朗日插值这边的题要么就是求个 \(S_k(n)\)，要么就是看出什么东西可以用多项式表示然后直接插出来系数接着做。
总的来说就是个小工具。
首先需要说明的是，\(f(x) - f(x - 1)\) 是一个关于 \(x\) 的 \(n\) 阶多项式，就表示 \(f(x)\) 是一个关于 \(x\) 的 \(n + 1\) 阶多项式。
简单挑几道题说说。
题面长，不粘。
“恰好”不好求，考虑容斥转化为“至少”。设 \(f_i\) 为至少 \(i\) 人被碾压的方案数。
钦定 \(i\) 个人被碾压，首先选出这些人来。然后对于每一科 \(j\) 我们在没有被钦定的 \(n - i -1\) 人里选择大于 B 神的 \(r_j - 1\) 个人，最后枚举 B 神在这一科的分数，然后每边随意选择。得到
发现后面那个对 \(k\) 的求和很难直接做，因此需要针对它化简一下。其实做差分后也能看出后面的是个对 \(u_i\) 的 \(n\) 次多项式。推一下：
好了。然后就没啥了，你直接上一个拉插就能在 \(O(n^2 m)\) 的复杂度内得到每个 \(f_i\)。答案即为
似乎有高论可以把对 \(k\) 的求和那部分优化到 \(O(\text M(a + b)\log n)\)。
令
则我们知道
答案就是对 \(F_a F_b\) 做远点求值。自然可以用 Bostan-Mori 方法做到 \(O(\text M(a + b) \log n)\)。将分母展开并用上指标反转可以只扫一遍 \(G_a G_b\) 的系数，组合数的维护好像要乘 \(O(\log \text{mod})\)。
冷静分析，求出 \(G_a G_b\) 是 \(O(\text{M}(a + b))\)，扫系数是 \(O((a + b) \log \text{mod})\)，快速幂的复杂度也是1log。还挺优？还要再乘 \(m\) 就是了。
总时间复杂度 \(O(nm \log u)\)。
code
胡写的。
#include <bits/stdc++.h>
using namespace std;
using pii = pair<int,int>; using vi = vector<int>; using vp = vector<pii>; using ll = long long;
using ull = unsigned long long; using db = double; using ld = long double; using lll = __int128_t;
#define multi int T; cin >> T; while ( T -- )
#define timer cerr << 1. * clock() / CLOCKS_PER_SEC << '\n';
#define iot ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define rep(i,s,t) for (register int i = (s), i##_ = (t) + 1; i < i##_; ++ i)
#define pre(i,s,t) for (register int i = (s), i##_ = (t) - 1; i > i##_; -- i)
#define eb emplace_back
#define pb pop_back
const int N = 1e6 + 10;
const int inf = 0x3f3f3f3f;
const ll infll = 0x3f3f3f3f3f3f3f3fll;
// int mod;
// const int mod = 10007;
// const int mod = 469762049, g = 3;
// const int mod = 998244353; // const int g = 3;
// const int mod = 1004535809, g = 3;
const int mod = 1e9 + 7;
// const int mod = 1e9 + 9;
// const int mod = 1e9 + 3579, bse = 131;
/* add / sub */ template<typename T1,typename T2>T1 add(T1 a,T2 b){return(a+=b)>=mod?a-mod:a;}template<typename T1,typename...Args>T1 add(T1 a,Args...b){return add(a,add(b...));}template<typename T1,typename T2>T1 sub(T1 a,T2 b){return(a-=b)<0?a+mod:a;}template<typename T1,typename...Args>T1 sub(T1 a,Args...b){return sub(a,add(b...));}template<typename T1,typename T2>void addi(T1&a,T2 b){(a+=b)>=mod?(a-=mod):true;}template<typename T1,typename...Args>void addi(T1&a,Args...b){addi(a,add(b...));}template<typename T1,typename T2>void subi(T1&a,T2 b){(a-=b)<0?(a+=mod):true;}template<typename T1,typename...Args>void subi(T1&a,Args...b){subi(a,add(b...));}
/* Fastmod / mul */ struct FastMod{int m;ll b;void init(int _m){m=_m;if(m==0)m=1;b=((lll)1<<64)/m;}FastMod(int _m){init(_m);}int operator()(ll a){ll q=((lll)a*b)>>64;a-=q*m;if(a>=m)a-=m;return a;}}Mod(mod);int mul(int a,int b){return Mod(1ll*a*b);}template<typename T1,typename T2>int mul(T1 a,T2 b){return Mod((long long)(1ll*a*b));}template<typename T,typename...Args>int mul(T a,Args...b){return mul(a,mul(b...));}template<typename T1,typename T2>void muli(T1&a,T2 b){a=Mod(1ll*a*b);}template<typename T1,typename...Args>void muli(T1&a,Args...b){muli(a,mul(b...));} // /* trivial multiple function(idk whether its useful*/ int mul(int a, int b) { return 1ll * a * b % mod; } template <typename T1, typename T2> int mul(T1 a, T2 b) { return (long long)(1ll * a * b) % mod; } template <typename T, typename ...Args> int mul(T a, Args ...b) { return mul(a, mul(b...)); }
/* qp fac C */ template<typename T1,typename T2>T1 qp(T1 a,T2 b){T1 ret=1;for(;b>0;a=1ll*a*a%mod,b>>=1)if(b&1)ret=mul(ret,a);return ret;}vi __fac({1,1}),__ifc({1,1}),__inv({0,1});inline void ___prep(int n){static int i=2;if(i<n)for(__fac.resize(n),__ifc.resize(n),__inv.resize(n);i<n;i++)__fac[i]=mul(i,__fac[i-1]),__inv[i]=mul((mod-mod/i),__inv[mod%i]),__ifc[i]=mul(__inv[i],__ifc[i-1]);}inline int fac(int x){return ___prep(x+1),__fac[x];}inline int ifc(int x){return ___prep(x+1),__ifc[x];}inline int inv(int x){return ___prep(x+1),__inv[x];}inline int C(int n,int m){if(n<m or n<0 or m<0)return 0;return mul(fac(n),ifc(m),ifc(n-m));}
/* sum of i^k */ int S(int n,int k){vector<int>__pre(k+4),__suf(k+4),__pw(k+4),__pri(k+4);vector<bool>__vis(k+4,0);__pw[1]=1;for(int i=2,cnt=0;i<=k+2;++i){if(!__vis[i])__pri[++cnt]=i,__pw[i]=qp(i,k);for(int j=1;j<=cnt and i*__pri[j]<=k+2;++j){__vis[i*__pri[j]]=1;__pw[i*__pri[j]]=mul(__pw[i],__pw[__pri[j]]);if(i%__pri[j]==0)break;}}rep(i,2,k+2)__pw[i]=add(__pw[i],__pw[i-1]);__pre[0]=__suf[k+3]=1;rep(i,1,k+2)__pre[i]=mul(__pre[i-1],(n-i+mod));pre(i,k+2,1)__suf[i]=mul(__suf[i+1],(n-i+mod));int tmp=0,ret=0;rep(i,1,k+2){if((k-i)&1)subi(ret,mul(__pw[i],__pre[i-1],__suf[i+1],ifc(i-1),ifc(k+2-i)));else addi(ret,mul(__pw[i],__pre[i-1],__suf[i+1],ifc(i-1),ifc(k+2-i)));}return ret;}
int qp(int a, int b, int p = mod) {
int ret = 1;
while (b) {
if (b & 1) ret = 1ll * ret * a % p;
a = 1ll * a * a % p;
b >>= 1;
} return ret;
}
using db = double;
struct cp {
db x, y;
cp(db real = 0, db imag = 0) : x(real), y(imag){};
cp operator+(cp b) const { return {x + b.x, y + b.y}; }
cp operator-(cp b) const { return {x - b.x, y - b.y}; }
cp operator*(cp b) const { return {x * b.x - y * b.y, x * b.y + y * b.x}; }
};
using vcp = vector<cp>;
namespace FFT {
const db pi = acos(-1);
vcp Omega(int L) {
vcp w(L); w[1] = 1;
for (register int i = 2; i < L; i <<= 1) {
auto w0 = w.begin() + i / 2, w1 = w.begin() + i;
cp wn(cos(pi / i), sin(pi / i));
for (register int j = 0; j < i; j += 2)
w1[j] = w0[j >> 1], w1[j + 1] = w1[j] * wn;
} return w;
} auto W = Omega(1 << 21);
void DIF(cp *a, int n) {
cp x, y;
for (register int k = n >> 1; k; k >>= 1)
for (register int i = 0; i < n; i += k << 1)
for (register int j = 0; j < k; ++j)
x = a[i + j], y = a[i + j + k],
a[i + j + k] = (x - y) * W[k + j], a[i + j] = x + y;
}
void IDIT(cp *a, int n) {
cp x, y;
for (register int k = 1; k < n; k <<= 1)
for (register int i = 0; i < n; i += k << 1)
for (register int j = 0; j < k; ++j)
x = a[i + j], y = a[i + j + k] * W[k + j],
a[i + j + k] = x - y, a[i + j] = x + y;
const db Inv = 1. / n;
rep(i, 0, n - 1) a[i].x *= Inv, a[i].y *= Inv;
reverse(a + 1, a + n);
}
}
namespace Polynomial {
using poly = vector<int>;
void DFT(vcp &a) { FFT::DIF(a.data(), a.size()); }
void IDFT(vcp &a) { FFT::IDIT(a.data(), a.size()); }
int norm(int n) { return 1 << (__lg(n - 1) + 1); }
poly operator*(poly &a, poly &b) {
int n = a.size() + b.size() - 1;
vcp c(norm(n));
rep(i, 0, a.size() - 1) c[i].x = a[i];
rep(i, 0, b.size() - 1) c[i].y = b[i];
DFT(c);
rep(i, 0, a.size() - 1) c[i] = c[i] * c[i];
IDFT(c), a.resize(n);
rep(i, 0, n - 1) a[i] = int(c[i].y * .5 + .5);
return a;
}
poly conv(const poly &a, const poly &b, const int&P = mod) {
int n = a.size(), m = b.size(), o = n + m - 1, l = norm(o);
vcp A(l), B(l), c0(l), c1(l);
for (register int i = 0; i < n; i++) A[i] = cp(a[i] & 0x7fff, a[i] >> 15);
for (register int i = 0; i < m; i++) B[i] = cp(b[i] & 0x7fff, b[i] >> 15);
FFT::DIF(A.data(), l), FFT::DIF(B.data(), l);
for (register int k = 1, i = 0, j; k < l; k <<= 1)
for ( ; i < k * 2; i++) {
j = i ^ k - 1;
c0[i] = cp(A[i].x + A[j].x, A[i].y - A[j].y) * B[i] * 0.5;
c1[i] = cp(A[i].y + A[j].y, -A[i].x + A[j].x) * B[i] * 0.5;
}
FFT::IDIT(c0.data(), l), FFT::IDIT(c1.data(), l);
poly res(o);
for (register int i = 0; i < o; i++) {
ll c00 = c0[i].x + 0.5, c01 = c0[i].y + 0.5, c10 = c1[i].x + 0.5, c11 = c1[i].y + 0.5;
res[i] = (c00 + ((c01 + c10) % P << 15) + (c11 % P << 30)) % P;
}
return res;
}
poly Inv(const poly& f, const int &p = mod) {
int o = f.size(), _n = norm(o); poly g(1), h;
g[0] = qp(f[0], p - 2, p);
for (int len = 2; len <= _n; len <<= 1) {
h.assign(f.begin(), f.begin() + min(len, (int)f.size()));
h = conv( g, conv( g, h, p ), p );
g.resize(len);
for (int i = 0; i < len; ++ i)
g[i] = (2ll * g[i] - h[i] + p) % p;
} g.resize(o);
return g;
}
poly Deri(const poly& f, const int& p = mod) {
poly ret; ret.resize(f.size());
for (int i = 0; i + 1 < f.size(); ++i) {
ret[i] = 1ll * f[i+1] * (i+1) % p;
} return ret;
} poly Intg(const poly& f, const int& p = mod) {
poly ret; ret.resize(f.size());
static int __Inv[N]; __Inv[0] = __Inv[1] = 1;
for (int i = 2; i < f.size(); ++ i)
__Inv[i] = 1ll * (p - p / i) * __Inv[p % i] % p;
for (int i = f.size() - 1; i >= 1; --i) {
ret[i] = 1ll * f[i-1] * __Inv[i] % p;
} return ret;
}
poly Ln(const poly& f, const int& p = mod) {
poly ret = Deri(f, p), Iv = Inv(f, p);
ret = conv(ret, Iv, p); ret = Intg(ret, p);
ret.resize(f.size());
return ret;
}
poly Exp(const poly& f, const int& p = mod) {
poly ret; ret.resize(1); ret[0] = 1;
poly A, B;
for (register int len = 2; len < (f.size() << 1); len <<= 1) {
ret.resize(len); A = Ln(ret);
B.assign(f.begin(), f.begin() + min(len, (int)f.size()));
if (B.size() < A.size()) B.resize(A.size());
for (int i = 0; i < A.size(); ++ i){
B[i] -= A[i]; if (B[i] < 0) B[i] += p;
} B[0]++; ret = conv(ret, B);
ret.resize(len);
} ret.resize(f.size());
return ret;
}
} using namespace Polynomial;
poly euler_l(int n) {
if (n == 0) {
poly A;
A.resize(1);
A[0] = 1;
return A;
}
poly A, B;
A.resize(n + 1);
B.resize(n + 1);
rep(i,0,n) A[i] = ((i & 1) ? Mod(mod - C(n + 1, i)) : C(n + 1, i));
rep(i,0,n) B[i] = qp(i + 1, n);
A = conv(A, B);
A.resize(n + 1);
pre(i,n,1) A[i] = A[i - 1];
A[0] = 0;
// cout << n << " : ";
// rep(i,1,n) cout << A[i] << ' ' ;
// cout << endl;
return A;
}
int G(int n, int a, int b) {
poly F = conv(euler_l(a), euler_l(b));
int ret = 0, now = 1;
rep(i,1,a+b+1) muli(now, n + 1 + i);
rep(i,0,a + b) {
now = mul(now, Mod(n - i + 1 + mod), qp(Mod(n - i + a + b + 2 + mod), mod - 2));
ret = add(ret, mul(now, F[i]));
// cout << i << ' ' << a + b << ' ' << F[i] << endl;
} now = 1;
rep(i,1,a+b+1) muli(now, i);
ret = mul(ret, qp(now, mod - 2));
return ret;
}
int n, m, k, a[105], s[105], ans;
int u[105], r[105];
signed main() {
// cin >> n;
// poly F = euler_l(n);
// for (auto v : F) cout << v << ' ';
// cout << endl;
cin >> n >> m >> k;
rep(i,1,m) cin >> u[i];
rep(i,1,m) cin >> r[i];
rep(j,1,m) {
/*
rep(t,0,r[j]-1) {
s[j] = add(s[j], mul( qp(u[j], t), C(r[j] - 1, t), ((r[j] - 1 - t) & 1) ? mod - 1 : 1, S(u[j], n - 1 - t) ));
}
*/
s[j] = G(u[j], n - r[j], r[j] - 1);
}
rep(i,k,n-1) {
a[i] = C(n - 1, i);
rep(j,1,m) {
a[i] = mul( a[i], C(n - i - 1, n - i - r[j]), s[j] );
}
}
rep(i,k,n-1)
ans = add(ans, mul(((i - k) & 1) ? mod - 1 : 1, C(i, k), a[i] ) );
cout << ans << '\n';
}
\(T\) 组询问，每组询问给定 \(k,a,n,d\)。
定义
\[f(n)=\sum_{i=1}^n i^k\qquad g(n)= \sum_{i=1}^n f(i) \qquad h(n) = \sum_{i=0}^n g(a + id) \]求
\[h(n) \bmod 1234567891 \]的值。
通过差分几次能得到 \(h(n)\) 是一个关于 \(n\) 的 \(k + 3\) 次多项式。
我们能方便的求得 \(f\)，自然可以方便地求得 \(g\) 的前 \(k + 2\) 个点值。自然可以方便地求得答案的前 \(k + 3\) 个点值。
然后一层一层地求即可。
一个长度为 \(n\)、元素值域在 \([1,k]\) 间且各个元素互不相同的序列 \(a\) 对答案的贡献为 \(\prod_{i=1}^n a_i\)。给定 \(k, n_0\) 对每个 \(1\le n \le n_0\)，求所有本质不同的长度为 \(n\) 的序列的权值和对 \(998244353\) 取模的值。
\(1\le n_0 \le 5\times 10^5,\ n_0\le k < 998244353\)。
谁用拉插做这种题啊，又难想复杂度也不优，还不如直接 DP 转生成函数！
我们考虑一个 DP 计算答案。元素互不相同启发我们在值域上做长度相关的 DP，设 \(f_{i, j}\) 为当前决策 \(i\) 选不选，当前计数的序列长度为 \(j\)。考虑当前值 \(i\) 是否选择，容易写出转移方程
然后对第一维求和得到
也就是
\(F_0(x) = 1\)。我们需要提取这玩意的第 \(n\) 项系数。容易想到一个 \(\prod\) 的经典套路是 \(\ln + \exp\)，首先作对数。
考虑到后面这个东西不是很好刻画，我们可以通过比对系数的方式来转化。具体地，我们考虑对后面的 \(\sum\) 单独构造生成函数，在需要转化回 \(\ln F_k(x)\) 的时候第 \(t\) 项乘入 \((-1)^{t + 1}/t\) 即可。
考虑到这个生成函数的第 \(i\) 项是等比数列求和的形式，我们可以通过 \(e^{px}\) 的 \(\text{EGF}\) 刻画其形式。具体地，我们构造
随后提取 \(x^t / t!\) 即可。带回原式得到
我们容易对给定的 \(t\) 构造出 \(e^{tx}\)，随后只需要做一次多项式求逆便可得到 \(G(x)\) 了。由于上下两个多项式的常数均为 \(0\)，因此可以作约分，最终多项式求逆是不受影响的。得到 \(G(x)\) 后转化为 \(\ln F_k(x)\) 是显然的。对结果做 \(\exp\) 后我们就构造出了 \(F_k(x)\)，其前 \(k\) 项即为答案。
总时间复杂度视实现方法而定，为 \(O(n\log n)/O(\frac{n\log^2 n}{\log \log n})\)。半在线卷积真香啊！
对于一个长度为 \(n\)、其中元素取值为 \([1,n]\) 内整数的序列 \(A = (A_1, A_2, \dots, A_n)\)，以及一个整数 \(i(1\le i \le n)\)，我们按如下方式定义一个长度为 \(10^{100}\) 的序列 \(B_i = (B_{i, 1}, B_{i, 2}, \dots, B_{i, 10^{100}})\)：
- \(B_{i, 1} = i\)。
- \(B_{i, j + 1} = A_{B_{i, j}} (1 \le j < 10^{100})\)。
此外，我们定义 \(S_i\) 为 \(B_i\) 中只出现一次的整数的数量。定义一个序列的贡献为 \(\sum_{i=1}^n S_i\)。
给定正整数 \(n, m\)，你需要求出总共 \(n^n\) 种可能的序列的贡献之和模 \(m\) 的值。
\(1\le n \le 2\times 10^5,\ 10^8 \le m\le 10^9\)。
由于答案具有对称性，任意 \(S_i\) 的贡献等于 \(S_1\)，接下来我们只需要讨论 \(S_1\)。下文中简记 \(B_{1, i}\) 为 \(B_i\)。
我们令 \(k\) 为一种可能中满足前 \(i\) 个元素彼此不同的最大的 \(i\)，随后来分别对确定的 \(k\) 计数种类数。
我们已经确定了 \(B_1 = 1\)，这样就需要从 \(n-1\) 个彼此不同的元素中选 \(k-1\) 个，也就有 \(A(n-1, k-1)\) 种可能性。关注 \(B_{k + 1}\)，可以发现其必定与前 \(k\) 个中的一个元素相同。假设 \(B_p = B_{k + 1} (1\le p \le k)\)，则我们可以发现，\(B_p\) 到 \(B_k\) 的子序列定会无限循环下去。
因此，只有 \(B_1\) 到 \(B_{p-1}\) 对应的子序列是只在 \(B\) 中出现一次的序列，对答案的贡献即为 \(p-1\)。对所有可能的子序列贡献求和得到这一部分对答案的贡献 \(\sum_{p=1}^k (p - 1) = p(p-1)/2\)。
现在固定 \(B_1\) 到 \(B_{k + 1}\) 段的序列，考虑有多少 \(A\) 满足当前的条件。由于 \(B_i\) 固定等价于 \(A_{B_{i-1}}\) 被确定，可以发现 \(A_{B_1}, A_{B_2}, \dots, A_{B_k}\) 都是已经确定了的，因此 \(A\) 中没有确定的元素还有 \(n - k\) 个。
因此对 \(k\) 的答案即为
对 \(1\le k\le n\) 求和即可得到 \(S_1\)。总时间复杂度 \(O(n)\)。
以下是博客签名，与正文无关。
请按如下方式引用此页：
本文作者 joke3579，原文链接：https://www.cnblogs.com/joke3579/p/chitchat230108.html。
遵循 CC BY-NC-SA 4.0 协议。
请读者尽量不要在评论区发布与博客内文完全无关的评论，视情况可能删除。
