牛客 11259 H Scholomance Academy 题解

传送门

数学板子大综合

数论+组合数学+多项式+线性递推板子套板子题

干脆加上群论和线性代数凑个整好了

---

【大意】

规定 (displaystyle G(N)=sum_{k_1+k_2+cdots+k_t=N}oldsymbol F(p_1^{k_1}p_2^{k_2}cdots p_t^{k_t}))

且 (displaystyle oldsymbol F(n)=sum_{a_1a_2cdots a_m}oldsymbol varphi(a_1)oldsymbol varphi(a_2)cdots oldsymbol varphi(a_m))

其中 (oldsymbol varphi(n)) 是 (n) 的欧拉函数值

现给定 (N, m, t, {p_1, p_2,cdots p_t})，求解 (G(N)mod 998244353)

---

【分析】

先从 (F) 下手，容易看出：

(displaystyle m=2 o oldsymbol F(n)=sum_{dmid n}oldsymbol varphi(d)oldsymbol varphi({nover d}))

数学归纳法归纳一下，发现 (oldsymbol F(n)=oldsymbol varphi^m(n))

此处的乘法为狄利克雷卷积

考虑到 (displaystyle oldsymbol varphi(p^k)=p^k(1-{1over p}), k>0) 且 (oldsymbol varphi(1)=oldsymbol varphi(p^0)=1)

所以有 (displaystyle oldsymbol varphi(p^k)=p^k(1-{1over p})^{[k>0]})

由于 (oldsymbol F=oldsymbol varphi^m) ，故 (oldsymbol F) 为积性函数，我们考虑其 (p^k) 的值

因此 (displaystyle oldsymbol F(p^k)=sum_{i_1+i_2+cdots +i_m=k}oldsymbol varphi(p^{i_1})oldsymbol varphi(p^{i_2})cdots oldsymbol varphi(p^{i_m}))

代入上式得 (displaystyle oldsymbol F(p^k)=p^k sum_{i_1+i_2+cdots +i_m=k}(1-{1over p})^{[i_1>0]+[i_2>0]+cdots [i_m>0]})

考虑枚举非零的个数 (r) ，剩下的 ((m-r)) 个全部为 (0) ，等价于将 (k) 个完全一样的小球放入 (r) 个完全一样的盒子，方案数为 (dbinom{k-1}{r-1})

代入得 (displaystyle oldsymbol F(p^k)=p^k sum_{r=1}^mdbinom m r dbinom{k-1}{r-1}(1-{1over p})^r,k>0)

---

现考虑 (oldsymbol F) 的积性，有：(displaystyle G(N)=sum_{k_1+k_2+cdots +k_m=N}oldsymbol F(p_1^{k_1})oldsymbol F(p_2^{k_2})cdots oldsymbol 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^nsum_{r=1}^mdbinom m r dbinom{n-1}{r-1}(1-{1over p})^r)

无法求和，故更换求和顺序，先求和 (r)

(displaystyle F_i(x)=1+sum_{r=1}^mdbinom m r (1-{1over p})^rsum_{n=1}^infty x^np_i^ndbinom{n-1}{r-1})

观察式子：(displaystyle sum_{n=1}^infty x^np_i^ndbinom{n-1}{r-1})

由于要 (ngeq r) 时 (dbinom {n-1}{r-1}) 才为正整数，且等于 (dbinom{n-1}{n-r})

故更改求和范围得：

(displaystyle sum_{n=1}^infty x^np_i^ndbinom{n-1}{r-1}=sum_{n=r}^infty x^np_i^n dbinom{n-1}{n-r}=x^rp_i^rsum_{n=r}^infty x^{n-r}p_i^{n-r}dbinom{(n-r)+r-1}{n-r})

将求和变量换元 (n-r o n) 得：

(displaystyle sum_{n=1}^infty x^np_i^ndbinom{n-1}{r-1}=x^rp_i^rsum_{n=0}^infty x^np_i^ndbinom{n+r-1}{n}=x^rp_i^rcdot {1over (1-xp_i)^r})

故将结果代回生成函数得：

(displaystyle F_i(x)=1+sum_{r=1}^mdbinom m r (1-{1over p_i})^rx^rp^rcdot {1over (1-xp_i)^r}=1+sum_{r=0}^m dbinom m r [(1-{1over p_i})cdot xcdot p_icdot {1over 1-xp_i}]^r)

发现是二项式定理，所以

(displaystyle F_i(x)=1+(1+{p_ix-xover 1-p_ix})^m-1=({1-xover 1-p_ix})^m)

所以代回生成函数，得到：

(displaystyle G(x)=prod_{i=1}^tF_i(x)=(prod_{i=1}^t {1-xover 1-p_ix})^m)

---

然而 (nleq 10^9) ，无法直接跑多项式求解。

而且由于其最短递推式长度 (k) 可能会达到 (10^5sim 10^6) 级别，甚至无法用 BM+(O(k^2log n)) 线性递推求解

因此这边需要用到其他方法求出最短递推式，而后用 (O(klog klog 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)={1over P(x)}=sum_{n=0}^infty r_nx^n)

由 (displaystyle 1=P(x)cdot R(x)) 得 (displaystyle forall n> m o 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 ngeq m+mu+1 o sum_{i=0}^m p_ig_{n-i}=sum_{i+j+k=n}p_iq_jr_k=sum_{j=0}^mu q_jsum_{i+k=n-j}p_ir_k)

由于 (n-jgeq 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_jsum_{i+k=n-j}p_ir_k=0) 在 (ngeq m+mu+1) 时恒成立

故移项得 (displaystyle g_n=-{1over p_0}sum_{i=1}^m p_ig_{n-i})

考虑到本题中 (displaystyle P(x)=prod_{i=1}^t(1-p_ix)^m o p_0=1) 以及线性递推的形式 (displaystyle a_n=sum_{i=1}^k f_ia_{n-i})

对应位置对应即可求解