牛客 11259 H Scholomance Academy 题解
数学板子大综合
数论+组合数学+多项式+线性递推板子套板子题
干脆加上群论和线性代数凑个整好了
【大意】
规定 \(\displaystyle G(N)=\sum_{k_1+k_2+\cdots+k_t=N}\boldsymbol F(p_1^{k_1}p_2^{k_2}\cdots p_t^{k_t})\)
且 \(\displaystyle \boldsymbol F(n)=\sum_{a_1a_2\cdots a_m}\boldsymbol \varphi(a_1)\boldsymbol \varphi(a_2)\cdots \boldsymbol \varphi(a_m)\)
其中 \(\boldsymbol \varphi(n)\) 是 \(n\) 的欧拉函数值
现给定 \(N, m, t, \{p_1, p_2,\cdots p_t\}\)，求解 \(G(N)\bmod 998244353\)
【分析】
先从 \(F\) 下手，容易看出：
\(\displaystyle m=2\to \boldsymbol F(n)=\sum_{d\mid n}\boldsymbol \varphi(d)\boldsymbol \varphi({n\over d})\)
数学归纳法归纳一下，发现 \(\boldsymbol F(n)=\boldsymbol \varphi^m(n)\)
此处的乘法为狄利克雷卷积
考虑到 \(\displaystyle \boldsymbol \varphi(p^k)=p^k(1-{1\over p}), k>0\) 且 \(\boldsymbol \varphi(1)=\boldsymbol \varphi(p^0)=1\)
所以有 \(\displaystyle \boldsymbol \varphi(p^k)=p^k(1-{1\over p})^{[k>0]}\)
由于 \(\boldsymbol F=\boldsymbol \varphi^m\) ，故 \(\boldsymbol F\) 为积性函数，我们考虑其 \(p^k\) 的值
因此 \(\displaystyle \boldsymbol F(p^k)=\sum_{i_1+i_2+\cdots +i_m=k}\boldsymbol \varphi(p^{i_1})\boldsymbol \varphi(p^{i_2})\cdots \boldsymbol \varphi(p^{i_m})\)
代入上式得 \(\displaystyle \boldsymbol F(p^k)=p^k \sum_{i_1+i_2+\cdots +i_m=k}(1-{1\over p})^{[i_1>0]+[i_2>0]+\cdots [i_m>0]}\)
考虑枚举非零的个数 \(r\) ，剩下的 \((m-r)\) 个全部为 \(0\) ，等价于将 \(k\) 个完全一样的小球放入 \(r\) 个完全一样的盒子，方案数为 \(\dbinom{k-1}{r-1}\)
代入得 \(\displaystyle \boldsymbol F(p^k)=p^k \sum_{r=1}^m\dbinom m r \dbinom{k-1}{r-1}(1-{1\over p})^r,k>0\)
现考虑 \(\boldsymbol F\) 的积性，有：\(\displaystyle G(N)=\sum_{k_1+k_2+\cdots +k_m=N}\boldsymbol F(p_1^{k_1})\boldsymbol F(p_2^{k_2})\cdots \boldsymbol F(p_t^{k_t})\)
形式上非常像生成函数
我们记 \(\displaystyle G(x)=\sum_{n=0}^\infty G(n)x^n, F_i(x)=\sum_{n=0}^\infty F(p_i^n)x^n\)
故 \(\displaystyle G=\prod_{i=1}^t F_i\)
因此我们只需要考虑求解 \(F_i(x)\)
而 \(\displaystyle F_i(x)=\sum_{n=0}^\infty F(p_i^n)x^n=1+\sum_{n=1}^\infty F(p_i^n) x^n\)
代入上一节的 \(F(p^k)\) 公式得
\(\displaystyle F_i(x)=1+\sum_{n=1}^\infty x^np_i^n\sum_{r=1}^m\dbinom m r \dbinom{n-1}{r-1}(1-{1\over p})^r\)
无法求和，故更换求和顺序，先求和 \(r\)
\(\displaystyle F_i(x)=1+\sum_{r=1}^m\dbinom m r (1-{1\over p})^r\sum_{n=1}^\infty x^np_i^n\dbinom{n-1}{r-1}\)
观察式子：\(\displaystyle \sum_{n=1}^\infty x^np_i^n\dbinom{n-1}{r-1}\)
由于要 \(n\geq r\) 时 \(\dbinom {n-1}{r-1}\) 才为正整数，且等于 \(\dbinom{n-1}{n-r}\)
故更改求和范围得：
\(\displaystyle \sum_{n=1}^\infty x^np_i^n\dbinom{n-1}{r-1}=\sum_{n=r}^\infty x^np_i^n \dbinom{n-1}{n-r}=x^rp_i^r\sum_{n=r}^\infty x^{n-r}p_i^{n-r}\dbinom{(n-r)+r-1}{n-r}\)
将求和变量换元 \(n-r\to n\) 得：
\(\displaystyle \sum_{n=1}^\infty x^np_i^n\dbinom{n-1}{r-1}=x^rp_i^r\sum_{n=0}^\infty x^np_i^n\dbinom{n+r-1}{n}=x^rp_i^r\cdot {1\over (1-xp_i)^r}\)
故将结果代回生成函数得：
\(\displaystyle F_i(x)=1+\sum_{r=1}^m\dbinom m r (1-{1\over p_i})^rx^rp^r\cdot {1\over (1-xp_i)^r}=1+\sum_{r=0}^m \dbinom m r [(1-{1\over p_i})\cdot x\cdot p_i\cdot {1\over 1-xp_i}]^r\)
发现是二项式定理，所以
\(\displaystyle F_i(x)=1+(1+{p_ix-x\over 1-p_ix})^m-1=({1-x\over 1-p_ix})^m\)
所以代回生成函数，得到：
\(\displaystyle G(x)=\prod_{i=1}^tF_i(x)=(\prod_{i=1}^t {1-x\over 1-p_ix})^m\)
然而 \(n\leq 10^9\) ，无法直接跑多项式求解。
而且由于其最短递推式长度 \(k\) 可能会达到 \(10^5\sim 10^6\) 级别，甚至无法用 BM+\(O(k^2\log n)\) 线性递推求解
因此这边需要用到其他方法求出最短递推式，而后用 \(O(k\log k\log n)\) 的方法求解线性递推
这边先给出结论：若 \(\displaystyle G(x)={Q(x)\over P(x)}\) ，其中 \(Q(x),P(x)\) 分别是 \(\mu,m\) 次多项式
则对于 \((m+\mu)\) 及以内的系数，需要暴力求解
而对于 \(deg=(m+\mu+1)\) 及以上的系数，由长度为 \(m\) 的线性递推求解
即 \((m+\mu)\) 及以内的系数无递推关系，以上的有长度为 \(m\) 的线性递推关系
该线性递推关系为 \(\displaystyle g_n=-\sum_{i=1}^m p_ig_{n-i}\)
证明见下文
故我们先使用分治 NTT 或多项式启发式合并，求解分母 \(P(x)\) ；用组合数展开求解分子 \(P(x)\)
然后对分母进行多项式取逆，暴力打出 \(G(x)\) 在模 \(x^{deg}\) 意义下的值
之后对于询问进行分支：
对于 \(deg\) 范围内的询问 \(n\) ，直接输出答案
对于 \(deg\) 范围外的询问 \(n\) ，丢进线性递推直接求解
这里强调一下，由于递推式从 \(deg\) 项开始有效，故递推关系 \(\displaystyle g_n=-\sum_{i=1}^m p_ig_{n-i}\) 的实际最小项为 \(g_{deg-n}\)
而线性递推的实际最小项为 \(g_0\)
因此需要对 \(n\) 减去偏移量 \((deg-n)\) ，计算线性递推后，乘上的向量也应是从 \(g_{deg-n}\) 开始的；即向量也需要加上一个 \((deg-n)\) 的偏移量
【代码】
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<ll, ll> pii;
typedef double db;
typedef vector<int> vi;
#define fi first
#define se second
#define rep(i, a, b) for(int i=(a); i<(b); ++i)
#define per(i, a, b) for(int i=(b)-1;i>=(a);--i)
#define pb push_back
#define mp make_pair
#define sz(a) (int)a.size()
#define all(a) a.begin(), a.end()
const int LimBit=19, M = 1<<LimBit<<1, P=998244353;
int a[M], b[M], c[M], d[M], moda[M], modb[M];
vi T[M];
int n, t, m, f[M], g[M], h[M];
inline int add(int a, int b) { return (a+=b)>=P?a-P:a; }
inline int dis(int a, int b) { return (a-=b)<0?a+P:a; }
inline int mul(int a, int b) { return (ll)a*b%P; }
ll kpow(ll a, int b){
ll c = 1;
for (; b; b >>= 1,a = a * a % P) if (b & 1) c = c * a %P;
return c;
}
struct NTT{
int N, na, nb, W[2][M+M], Rev[M+M], *w[2], *rev, invN, InV[M];
NTT(){
w[0]=W[0], w[1]=W[1], rev=Rev;
w[0][0]=w[1][0]=1;
for(int i=1,x=kpow(3,P>>LimBit),y=kpow(x,P-2); !(i>>LimBit); ++i){
rev[i]=(rev[i>>1]>>1)|((i&1)<<LimBit-1);
w[0][i]=(ll)w[0][i-1]*x%P, w[1][i]=(ll)w[1][i-1]*y%P;
}
int *lstw[2]={w[0], w[1]};
w[0]+=1<<LimBit, w[1]+=1<<LimBit, rev+=1<<LimBit;
for(int d=LimBit-1;d>=0;--d){
for(int i=0;!(i>>d);++i){
rev[i]=(rev[i>>1]>>1)|((i&1)<<d-1);
w[0][i]=lstw[0][i<<1], w[1][i]=lstw[1][i<<1];
}
lstw[0]=w[0], lstw[1]=w[1];
w[0]+=1<<d, w[1]+=1<<d, rev+=1<<d;
}
}
inline void work(){
w[0]=W[0]; w[1]=W[1]; rev=Rev;
for(int d=LimBit;1<<d>N;--d)
w[0]+=1<<d, w[1]+=1<<d, rev+=1<<d;
invN=kpow(N,P-2);
}
inline void FFT(int *a,int f){
for(int i=0;i<N;++i) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int i=1;i<N;i<<=1)
for(int j=0,t=N/(i<<1);j<N;j+=i<<1)
for(int k=0,l=0,x,y;k<i;k++,l+=t)
x=(ll)w[f][l]*a[j+k+i]%P, a[j+k+i]=dis(a[j+k],x), a[j+k]=add(a[j+k],x);
if(f) for(int i=0;i<N;++i) a[i]=(ll)a[i]*invN%P;
}
inline void doit(int *a,int *b,int na,int nb){//3*NTT
for(N=1;N<na+nb-1;N<<=1);
for(int i=na;i<N;++i) a[i]=0;
for(int i=nb;i<N;++i) b[i]=0;
work(); FFT(a,0); FFT(b,0);
for(int i=0;i<N;++i) a[i]=(ll)a[i]*b[i]%P;
FFT(a,1);
}
inline void get_inv(int *f,int *g,int n){//6*NTT
g[0]=kpow(f[0],P-2);
int l;
for(l=1;l<n;l<<=1){
for(int i=0;i<l;++i) a[i]=f[i];
for(int i=l;i<l<<1;++i) a[i]=f[i],g[i]=0;
N=l<<2;
for(int i=l<<1;i<N;++i) a[i]=g[i]=0;
work(); FFT(a,0); FFT(g,0);
for(int i=0;i<N;++i) g[i]=(ll)g[i]*dis(2,(ll)g[i]*a[i]%P)%P;
FFT(g,1);
}
for(int i=n;i<N;++i) g[i]=0;
}
priority_queue<pii> H;
inline void Merge(int cnt){
while(H.size()) H.pop();
for(int i=1;i<=cnt;++i) H.emplace( pii(-T[i].size(), i) );
while(H.size()>=2){
int x=H.top().se; H.pop();
int y=H.top().se; H.pop();
for(int i=0;i<T[x].size();++i) a[i]=T[x][i];
for(int i=0;i<T[y].size();++i) b[i]=T[y][i];
doit(a, b, T[x].size(), T[y].size());
int na=T[x].size()+T[y].size();
T[x].clear(); T[y].clear();
for(int i=0;i<na;++i) T[x].push_back(a[i]);
H.emplace( pii(-T[x].size(), x) );
}
swap(T[H.top().se], T[1]);
H.pop();
}
inline void get_mod(int *f,int *g,int n,int m){
while(n>=0&&!f[n-1]) --n;
if(n<m) return ;
reverse(g, g+m);
get_inv(g, modb, n-m+1);
reverse(g, g+m);
for(int i=0;i<n;++i) moda[i]=f[i];
reverse(moda, moda+n);
doit(moda, modb, n-m+1, n-m+1);
reverse(moda, moda+n-m+1);
for(int i=0;i<m;++i) modb[i]=g[i];
doit(moda, modb, n-m+1, m);
for(int i=0;i<n;++i) f[i]=dis(f[i], moda[i]);
}
inline void get_xpow(int *g,int *m,int x,int nm){
if(x==0){
g[0]=1;
return ;
}
get_xpow(g, m, x>>1, nm);
for(N=1;N<nm+nm-3;N<<=1);
for(int i=nm-1;i<N;++i) g[i]=0;
work(); FFT(g, 0);
for(int i=0;i<N;++i) g[i]=(ll)g[i]*g[i]%P;
FFT(g, 1);
if(x&1){
for(int i=N;i>=1;--i) g[i]=g[i-1];
g[0]=0;
get_mod(g, m, N+1, nm);
}
else get_mod(g, m, N, nm);
for(N=1;N<nm+nm-3;N<<=1);
for(int i=nm-1;i<N;++i) g[i]=0;
}
} ntt;
const int MAXN=1e6+1;
int Inv[MAXN], deg;
inline void init(){
cin>>n>>t>>m;
for(int i=1, p;i<=t;++i){
cin>>p;
T[i].push_back(1);
T[i].push_back( dis(0, p) );
}
ntt.Merge(t);
for(int i=2;i<=m;++i) T[i]=T[1];
ntt.Merge(m);
deg=T[1].size()+m*t;
for(int i=0;i<deg;++i) g[i]=(i<T[1].size()?T[1][i]:0);
ntt.get_inv(g, f, deg);
Inv[1]=1;
for(int i=2;i<MAXN;++i)
Inv[i]=P-(ll)Inv[P%i]*(P/i)%P, g[i]=0;
for(int i=0, sgn=0, comb=1, I=m*t;i<=I;++i, sgn^=1){
if(sgn) g[i]=dis(0, comb);
else g[i]=comb;
comb=(ll)comb*(I-i)%P*Inv[i+1]%P;
}
ntt.doit(f, g, deg, m*t+1);
}
struct LinearRecurrence{
int k, g[M], c[M];
inline void pre(const vector<int> &f) {
k=f.size()-1;
for(int i=0;i<k;++i) g[i]=f[k-i];
g[k]=1;
}
inline int ans(int *a, int n) {
ntt.get_xpow(c, g, n, k+1);
int res=0;
for(int i=0;i<k;++i)
res=add(res, (ll)c[i]*a[i]%P);
return res;
}
inline int work(int *a, int n){
int offset=deg-k;
if(n<deg) return a[n];
else return ans(a+offset, n-offset);
}
}LR;
int main(){
ios::sync_with_stdio(0);
cin.tie(0); cout.tie(0);
init();
LR.pre(T[1]);
cout<<LR.work(f, n);
cout.flush();
return 0;
}
【证明】
感谢 Stump 贡献的方向
设 \(\displaystyle P(x)=\sum_{n=0}^m p_nx^n, Q(x)=\sum_{n=0}^\mu q_nx^n\)
设 \(\displaystyle R(x)={1\over P(x)}=\sum_{n=0}^\infty r_nx^n\)
由 \(\displaystyle 1=P(x)\cdot R(x)\) 得 \(\displaystyle \forall n> m\to \sum_{i=0}^m p_ir_{n-i}=\sum_{i+j=n}p_ir_j=0\)
由 \(\displaystyle G(x)={Q(x)\over P(x)}=Q(x)\cdot R(x)\) 得 \(\displaystyle g_n=\sum_{i+j=n}q_ir_j\)
故 \(\displaystyle \forall n\geq m+\mu+1\to \sum_{i=0}^m p_ig_{n-i}=\sum_{i+j+k=n}p_iq_jr_k=\sum_{j=0}^\mu q_j\sum_{i+k=n-j}p_ir_k\)
由于 \(n-j\geq n-\mu \geq m+1>m\) 故 \(\displaystyle \sum_{i+k=n-j}p_ir_k=0\) 恒成立
因此 \(\displaystyle \sum_{i=0}^m p_ig_{n-i}=\sum_{j=0}^\mu q_j\sum_{i+k=n-j}p_ir_k=0\) 在 \(n\geq m+\mu+1\) 时恒成立
故移项得 \(\displaystyle g_n=-{1\over p_0}\sum_{i=1}^m p_ig_{n-i}\)
考虑到本题中 \(\displaystyle P(x)=\prod_{i=1}^t(1-p_ix)^m\to p_0=1\) 以及线性递推的形式 \(\displaystyle a_n=\sum_{i=1}^k f_ia_{n-i}\)
对应位置对应即可求解
