题解 P5401 【[CTS2019]珍珠】

djq 提到了一种时间复杂度更优的解法，看到洛谷题解区还没有，就来复读一遍吧（

省略亿点步骤，题目中要求的可以转化为如下问题：

\frac{n!}{2^d}[x^n]\sum_{i=0}^k \binom di(\text e^x-\text e^{-x})^i(\text e^x+\text e^{-x})^{d-i}

其中 k = \min(d,n-2m)。

由于 EGF 不太好搞，而关于 x 的式子只有 \text e^x 这样的，于是考虑转为 OGF：

F(x)=\sum_{i=0}^k \binom di(x-x^{-1})^i(x+x^{-1})^{d-i} =x^{-d}\sum_{i=0}^k \binom di (x^2-1)^i(x^2+1)^{d-i}

这样要求的就是 [x^n]F(\text e^x)，但形如 (x^2+1) 的形式还是不好看，改为：

G(x)=\sum_{i=0}^k \binom di(x-1)^i(x+1)^{d-i}

然后答案就能写成

\frac{n!}{2^d}[x^n]G(\text e^{2x})\text e^{-dx}

只要求出 G(x) 的系数，那么其卷积的 n 次项就能直接求出。
具体来讲，答案可以写为

2^{-d}\sum_{i=0}^dg_i (2i-d)^n

线性筛求幂，是 \Theta(d \log_d n) 的。

然后使用万能的（并不）求导大法，就能推出：

G'(x)=\sum_{i=0}^k \binom di(i(x-1)^{i-1}(x+1)^{d-i}+(d-i)(x-1)^i(x+1)^{d-i-1}) 2dG(x)-2xG'(x)=d\binom{d-1}{k}((x-1)^k(x+1)^{d-k}-(x-1)^{k+1}(x+1)^{d-k-1}) dG(x)-xG'(x)=d\binom{d-1}{k}(x-1)^k(x+1)^{d-k-1}

设 H(x)=dG(x)-xG'(x)，就可以直接提取系数 (d-n)g_n=h_n 计算 G(x) 的系数（注意，h_d 为零，而 g_d 显然不为零，所以要单独算一下）。

要算 H(x) 很简单，只要能快速求形如

Q(x)=(x+1)^a(x-1)^b

的系数就好了。

Q'(x)=a(x+1)^{a-1}(x-1)^b+b(x+1)^a(x-1)^{b-1} = \frac{aQ(x)}{x+1}+\frac{bQ(x)}{x-1} (x^2-1)Q'(x)=a(x-1)Q(x)+b(x+1)Q(x) -(n+1)q_{n+1}=a(q_{n-1}-q_n)+b(q_{n-1}+q_n)-(n-1)q_{n-1}

显然可以 \Theta(d) 递推计算前 d 项。

总时间复杂度 \Theta(d \log_d n)。

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#define N 100003
#define ll long long
#define reg register
#define p 998244353
using namespace std;

inline int add(const int& x,const int& y){ return x+y>=p?x+y-p:x+y; }

inline int power(int a,int t){
    int res = 1;
    while(t){
        if(t&1) res = (ll)res*a%p;
        a = (ll)a*a%p;
        t >>= 1;
    }
    return res;
}

int inv[N],pw[N],pr[N>>2],fac[N],ifac[N];

void sieve(int n,int k){
    int cnt = 0;
    fac[0] = fac[1] = ifac[0] = ifac[1] = inv[1] = pw[1] = 1;
    for(reg int i=2;i<=n;++i){
        inv[i] = (ll)(p-p/i)*inv[p%i]%p;
        fac[i] = (ll)fac[i-1]*i%p;
        if(!pw[i]){
            pr[++cnt] = i;
            pw[i] = power(i,k);
        }
        for(reg int j=1;j<=cnt&&i*pr[j]<=n;++j){
            pw[i*pr[j]] = (ll)pw[i]*pw[pr[j]]%p;
            if(i%pr[j]==0) break;
        }
    }
    ifac[n] = power(fac[n],p-2);
    for(reg int i=n-1;i;--i) ifac[i] = (ll)ifac[i+1]*(i+1)%p;
}

inline int binom(int n,int m){ return (ll)fac[n]*ifac[m]%p*ifac[n-m]%p; }

int d,n,k,ans,sgn,a,b,bm;
int g[N],q[N];

int main(){
    scanf("%d%d%d",&d,&n,&k);
    if((k<<1)<=n-d){
        printf("%d",power(d,n));
        return 0;
    }else if((k<<1)>n){
        putchar('0');
        return 0;
    }
    k = min(d,n-(k<<1));
    sgn = n&1?-1:1;
    sieve(d+1,n);
    bm = (ll)d*binom(d-1,k)%p;
    a = d-k-1,b = k;
    q[0] = b&1?-1:1,q[1] = (b&1?-1:1)*a+((b-1)&1?-b:b);
    for(reg int i=1;i!=d;++i) q[i+1] = -((ll)a*(q[i-1]-q[i])%p+(ll)b*(q[i-1]+q[i])-(ll)(i-1)*q[i-1])%p*inv[i+1]%p;
    for(reg int i=0;i!=d;++i) g[i] = (ll)q[i]*bm%p*inv[d-i]%p;
    for(reg int i=0;i<=k;++i) g[d] = add(g[d],binom(d,i));
    for(reg int i=0;i<=d;++i){
        if((i<<1)<d) ans = (ans+(ll)g[i]*sgn*pw[d-(i<<1)])%p;
        else ans = (ans+(ll)g[i]*pw[(i<<1)-d])%p;
    }
    ans = (ll)ans*power(2,p-1-d)%p;
    printf("%d",ans<0?ans+p:ans);
    return 0;   
}