HDU 6900 Residual Polynomial【分治 NTT】

 

HDU 6900 Residual Polynomial 

题意：

给出一个多项式(f_1(x) = sum_{i=0}^na_ix^i)

对于任意(i>=2)，满足(f_i(x) = b_i(f_{i-1}(x))'+c_if_{i-1}(x))

要求得到(f_n(x))的各次项系数模(998244353)

(nle 10^5, 0leq a_i,b_i,c_i < 998244353)

题解：

考虑把(f_1,f_2,cdots,f_n)写成(n)列，其中(f_{ij})表示(f_i)的(j)次项的系数：

[egin{array}{cccc} f_{1,0} & f_{2,0} & f_{3,0} & cdots & f_{n0,} \ f_{1,1} & f_{2,1} & f_{3,1} & cdots & f_{n,1}\ vdots & vdots & vdots &ddots & vdots \ f_{1,n} & f_{2,n} & f_{3,n} & cdots & f_{n,n} end{array} ]

其中(f_{1,i}=a_i)

考虑类似(dp)的状态转移，那么可以发现， 对于(f_{ij})存在两种转移：

[egin{cases} f_{i,j}stackrel{c_{i+1}}{longrightarrow}f_{i+1,j} & i<n \ f_{i,j}stackrel{jcdot b_{i+1}}{longrightarrow} f_{i+1,j-1} & i<n , j>0 end{cases} ]

其中箭头表示乘自身然后加到右边

那么如果把转移看作边，可以发现每个状态（除了边界）向右连了一条边，向右上连了一条边

我们来考虑(f_{1,i})对(f_{n,j})的贡献，其中(ige j)，可以发现(dp)的转移其实就是一条条从(f_{1,i})到(f_{n,j})的路径，那么显然可以把所有路径单独分开来看，那么(f_{1,i})对(f_{n,j})的贡献为：(f_{1,i}cdot sum (prod operatorname{pathvalue}))，也就是从(f_{1,i})到(f_{n,j})的所有可行路径的边权乘积的和

先不考虑第二类转移中(jcdot b_{i+1})的(j)，那么对于任何转移(f_{1,i} ightarrow f_{n,j})，其实就是对于每个(2le kle n)，选择(b_k)或者(c_k)，其中(b_k)选择(i-j)个，(c_k)选择(n-1-(i-j))个，乘起来然后再把所有方案加起来

我们令选(x)个(b_k)和(n-1-x)个(c_k)的所有方案的和为(F(x))，那么(F(x))是可以用分治+(FFT)来得到

令(F(l,r,x))表示在区间([l,r))中选(x)个(b_k)的所有方案的和，那么可以得到(F(l,r,x)=sum_{i+j=x}F(l,mid,i)cdot F(mid,r,j))，其中(mid = lfloor frac{l+r}2 floor)

再算上之前没考虑的(j)的贡献乘积，可以得到(f_{n,j}=sum_{i-k=j}F(k)cdot f_{1,i}cdot frac{(i-1)!}{(j-1)!})

考虑把数组反向，也就是(f_{i,j})和(f_{i,n-j})互换，那么就可以得到(f_{n,j} = sum_{i+k=j}F(k)cdot f_{1,i}cdot frac{j!}{i!} = j!sum_{i+k=j}F(k)cdot frac{f_{1,i}}{i!})

那么就可以再做一次(FFT)就能得到答案了

view code

//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
const int MOD = 998244353;
const int FFTN = 1<<19;
const int MAXN = 2e5+7;
#define poly vector<int>
typedef unsigned long long int ull;
int ksm(int a, int b){
    int ret = 1;
    while(b){
        if(b&1) ret = 1ll * ret * a % MOD;
        b >>= 1;
        a = 1ll * a * a % MOD;
    }
    return ret;
}
namespace FFT{
    int w[FFTN+5],W[FFTN+5],R[FFTN+5];
    void FFTinit(){
        W[0]=1;
        W[1]=ksm(3,(MOD-1)/FFTN);
        for(int i = 2; i <= FFTN; i++) W[i]=1ll*W[i-1]*W[1]%MOD;
    }
    int FFTinit(int n){
        int L=1;
        for (;L<=n;L<<=1);
        for(int i = 0; i <= L - 1; i++) R[i]=(R[i>>1]>>1)|((i&1)?(L>>1):0);
        return L;
    }
    int A[FFTN+5],B[FFTN+5];
    ull p[FFTN+5];
    void DFT(int *a,int n){
        for(int i = 0; i < n; i++) p[R[i]]=a[i];
        for(int d = 1; d < n; d <<= 1){
            int len=FFTN/(d<<1);
            for(int i = 0, j = 0; i < d; i++, j += len) w[i]=W[j];
            for(int i = 0; i < n; i += (d<<1))	
                for (int j = 0; j < d; j++){
                    int y=p[i+j+d]*w[j]%MOD;
                    p[i+j+d]=p[i+j]+MOD-y;
                    p[i+j]+=y;
                }
            if (d==1<<15)
                for(int i = 0; i < n; i++) p[i]%=MOD;
        }
        for(int i = 0; i < n; i++) a[i]=p[i]%MOD;
    }
    void IDFT(int *a,int n){
        for(int i = 0; i < n; i++) p[R[i]]=a[i];
        for (int d=1;d<n;d<<=1){
            int len=FFTN/(d<<1);
            for (int i=0,j=FFTN;i<d;i++,j-=len) w[i]=W[j];
            for (int i=0;i<n;i+=(d<<1))	
                for (int j=0;j<d;j++){
                    int y=p[i+j+d]*w[j]%MOD;
                    p[i+j+d]=p[i+j]+MOD-y;
                    p[i+j]+=y;
                }
            if (d==1<<15)
                for(int i = 0; i < n; i++) p[i]%=MOD;
        }
        int val=ksm(n,MOD-2);
        for(int i = 0; i < n; i++) a[i]=p[i]*val%MOD;
    }
    poly Mul(const poly &a,const poly &b){
        int sza=a.size()-1,szb=b.size()-1;
        poly ans(sza+szb+1);
        if (sza<=30||szb<=30){
            for(int i = 0; i <= sza; i++) for(int j = 0; j <= szb; j++)
                ans[i+j]=(ans[i+j]+1ll*a[i]*b[j])%MOD;
            return ans; 
        }
        int L=FFTinit(sza+szb);
        for(int i = 0; i < L; i++) A[i]=(i<=sza?a[i]:0);
        for(int i = 0; i < L; i++) B[i]=(i<=szb?b[i]:0);
        DFT(A,L); DFT(B,L);
        for(int i = 0; i < L; i++) A[i]=1ll*A[i]*B[i]%MOD;
        IDFT(A,L);
        for(int i = 0; i <= sza + szb; i++) ans[i]=A[i];
        return ans; 
    }
}
int fac[MAXN], rfac[MAXN], inv[MAXN];
poly divide(int l, int r, vector<int> &B, vector<int> &C){ return l + 1 == r ? poly({C[l],B[l]}) : FFT::Mul(divide(l,(l+r)>>1,B,C), divide((l+r)>>1,r,B,C)); }
void solve(){
    int n; scanf("%d",&n);
    vector<int> A(n+1), B(n-1), C(n-1);
    for(int &x : A) scanf("%d",&x);
    for(int &x : B) scanf("%d",&x);
    for(int &x : C) scanf("%d",&x);
    poly f = divide(0,n-1,B,C);
    reverse(A.begin(),A.end());
    for(int i = 0; i <= n; i++) A[i] = 1ll * A[i] * fac[n-i] % MOD;
    poly W = FFT::Mul(A,f);
    for(int i = 0; i <= n; i++) W[i] = 1ll * W[i] * rfac[n-i] % MOD;
    for(int i = n; i >= 0; i--) printf("%d%c",W[i]," 
"[!i]);
}
int main(){
    fac[0] = rfac[0] = inv[1] = 1;
    for(int i = 1; i < MAXN; i++) fac[i] = 1ll * fac[i-1] * i % MOD;
    for(int i = 2; i < MAXN; i++) inv[i] = 1ll * (MOD - MOD / i) * inv[MOD % i] % MOD;
    for(int i = 1; i < MAXN; i++) rfac[i] = 1ll * rfac[i-1] * inv[i] % MOD;
    FFT::FFTinit();
    int tt; for(scanf("%d",&tt); tt--; solve());
    return 0;
}