题解 P4705 【玩游戏】

· · 题解

基础前置技能

1.基础数学
2.NTT(FFT)
3.多项式求逆
4*.Taylor展开

设 ans_k，为

ans_k=\sum_{i=1}^n \sum_{j=1}^m (a_i+b_j)^k

设 Ans(x) 为 ans 的生成函数

答案即为求多项式 Ans(x) * \frac{1}{nm}

然后....首先这个式子....

ans_k= \sum_{i=1}^n \sum_{j=1}^m (a_i+b_j)^k

考虑二项式展开

ans_k= \sum_{i=1}^n \sum_{j=1}^m \sum_{t=0}^k C_k^t a_i^{k-t} b_j^t

考虑交换求和顺序，把 k 提到前面

ans_k= \sum_{t=0}^k \sum_{i=1}^n \sum_{j=1}^m C_k^t a_i^{k-t} b_j^t ans_k= \sum_{t=0}^k C_k^t \sum_{i=1}^n a_i^{k-t}\sum_{j=1}^m b_j^t

可以发现，如果计两个新的多项式为 A(x),B(x)，令

A(x)=\sum_{i=1}^n a_i^x B(x)=\sum_{i=1}^m b_i^x

就可以得到一个很优美的化简:

ans_k= \sum_{t=0}^k C_k^t\ A(k-t)\ B(t)

组合数直接拆出来

ans_k= \sum_{t=0}^k \frac{k!}{t!(k-t)!}\ A(k-t)\ B(t) =k! \sum_{t=0}^k \frac{A(k-t)}{(k-t)!}\ \frac{B(t)}{t!}

非常好。现在的问题变成了A(x)=\sum_{i=1}^n a_i^x这个看起来非常不友好的多项式怎么求。(说实话不会的话真的是没法想出来的...)(Drench说是非常经典的问题)

先考虑一个多项式，即为F(x)，为

F(x)=\prod_{i=1}^n(a_ix+1)

这个东西可以分治求出来。具体一点大概就是...

F_l(x)=\prod_{i=1}^{(n/2)}\ (a_ix+1) F_r(x)=\prod_{i=(n/2)+1}^{n}(a_ix+1) F(x)=F_l(x)F_r(x)

求出这个东西之后，再对这个东西取ln....

ln(F(x))=G(x)

这是个多项式求ln，具体来说...

(ln(x))'=\frac{1}{x}

链式法则：如果

h(x)=f(g(x))

那么就会有

h'(x)=f'(g(x))g'(x)

这里的ln(x)就相当于f(x)

因为G(x)=f(F(x))

G'(x)=f'(F(x))F'(x)=\frac{F'(x)}{F(x)}

可能需要反应一下

然后，

G(x)=\int \frac{F'(x)}{F(x)}

至于多项式的求导和积分，其实非常简单，求导的时候a_{i}=(i+1)a^{(i+1)}, 积分的时候a_i=inv(i)a_{i-1}。关键是求出\frac{1}{F(x)}是一个多项式求逆的过程。P4238 【模板】多项式求逆

这么费劲求出 G(x) 干嘛呢？

......

泰勒(taylor)展开：

首先约定用f^{(n)}(x)表示f(x)的n阶导数。

若函数 f(x) 在包含 x0 的某个闭区间 [a,b] 上具有 n 阶导数，且在开区间 (a,b) 上具有 (n+1) 阶导数，则对闭区间 [a,b] 上任意一点 x ，成立下式：

上面是百度百科抄的(顺手引用了图片。)。

用非常不准确的人话来说，就是有一个函数f(x)，这个函数随便是什么log,sin,cos，只要这个f(x)可在任意点导任意次，我们就可以用一个非常奇怪的无穷项的多项式g(x)从某一个点x_0来拟合这样的一个函数。g(x)的每一项a_i为：

\frac{f^{(n)}(x_0)}{n!}

合起来就是

f(x)=\sum_{k=1}^{+\infty} \frac{f^{(k)}(x_0)}{k!}(x-x_0)^k

这个东西的大致思路就是让g(x_0)=f(x_0)，让g'(x_0)=f'(x_0)，让g''(x_0)=f''(x_0)，直到让g^{(+\infty)}(x_0)=f^{(+\infty)}(x_0)，这两个函数就可以说一样了。满足这个条件的g(x)写出来就是泰勒展开。当然没有这么简单，关于taylor展开知乎上有讲的非常好的，如果想学建议自己去看一看。

回到这个东西 F(x)=\prod_{i=1}^n(a_ix+1) 。用分治NTT计算出 F(x) 之后再对它多项式求ln得到 G(x) :

G(x)=ln(\ \ \prod_{i=1}^n(a_ix+1)\ \ )

非常明显的化简

G(x)=\sum_{i=1}^nln(a_ix+1)

然后我们考虑把里面的这个函数 ln(x) 在 1 这个点进行Taylor展开。

原始的公式

f(x)=\sum_{k=1}^{+\infty} \frac{f^{(k)}(x_0)}{k!}(x-x_0)^k $$(ln(x))'=x^{-1}\ \ \ \ \ \ \ (x^a)'=(a-1)x^{a-1}$$ $$(ln(x))^{(n)}=(-1)^{(n+1)}n!x^{-n}$$ 带入 $$ln(a_ix+1)=\sum_{k=1}^{+\infty} \frac{(-1)^{(k+1)}(k-1)!}{k!}(a_ix)^k$$ $$ln(a_ix+1)=\sum_{k=1}^{+\infty} \frac{(-1)^{(k+1)}}{k}a_i^kx^k$$ 因为在求$G(x)$的时候用到了多项式求逆，对$x^{K+1}$取模后，这个式子中大于$K$次的项没有意义. 对于单独的一项 $$ln(a_ix+1)=\sum_{k=1}^{K} \frac{(-1)^{(k+1)}}{k}a_i^kx^k$$ 所以 $$G(x)=\sum_{i=1}^nln(a_ix+1)=\sum_{i=1}^n\sum_{k=1}^{K} \frac{(-1)^{(k+1)}}{k}a_i^kx^k$$ 为了表达的清楚一点随便找一个$t$把$t(k)$记为$\frac{(-1)^{(k+1)}}{k} G(x)=\sum_{i=1}^nln(a_ix+1)=\sum_{i=1}^n\sum_{k=1}^{K} t(k)a_i^kx^k

显然可以交换一下顺序

G(x)=\sum_{k=1}^{K} t(k)\sum_{i=1}^na_i^kx^k G(x)=\sum_{k=1}^{K} t(k)A(k)x^k

哇！

然后我的代码写的巨型丑，常数莫名巨型大，卡了一次才A,仅供对拍...

#include <bits/stdc++.h>
#define rap(i,s,n) for(int i=s;i<=n;i++)
#define drap(i,n,s) for(int i=n;i>=s;i--)
#define N 410000
#define P 998244353
#define lb(x) ((x&(-x)))
#define ll long long
#define m(s,k) memset(s,k,sizeof s)
using namespace std;
char xB[1<<15],*xS=xB,*xTT=xB;
#define getc() (xS==xTT&&(xTT=(xS=xB)+fread(xB,1,1<<15,stdin),xS==xTT)?0:*xS++)
#define isd(c) ((c>='0'&&c<='9')||(c=='-'))
template<typename T>
inline bool rd(T& xa){
    char xchh; T f=1; while(xchh=getc(),(!isd(xchh))&&(xchh!=0));
    if(xchh==0) return 0; if(xchh=='-') xchh=getc(),f=-1; xa=xchh-'0';
    while(xchh=getc(),isd(xchh)) xa=xa*10+xchh-'0'; xa*=f; return 1;
}
ll mpow(ll a,ll k,ll p){ll res=1; while(k){if(k&1) res=res*a%p; k>>=1; a=a*a%p;} return res%p;}
ll Inv(ll x){return mpow(x,P-2,P);}
ll fac[N],inv[N],finv[N];
void Finit(int n){
    fac[0]=1; rap(i,1,n) fac[i]=fac[i-1]*i%P;
    inv[1]=1; rap(i,2,n) inv[i]=inv[P%i]*(P-P/i)%P;
    finv[0]=1; rap(i,1,n) finv[i]=finv[i-1]*inv[i]%P;
}
namespace Poly{
    ll c[N],f[N],f1[N],fi[N];
    int rev[N],lim,l;
    void init(int n){
        lim=1,l=0; while(lim<=n) lim<<=1,l++;
        rap(i,1,lim-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
    }
    void NTT(ll *A,int k){
        rap(i,1,lim-1) if(i<rev[i]) swap(A[i],A[rev[i]]); for(int l=1;l<lim;l<<=1){
            ll gn=mpow(3,(P-1)+k*(P-1)/l/2,P); for(int j=0;j<lim;j+=l*2){
                ll g=1; for(int k=0;k<l;k++,g=g*gn%P){
                    ll x=A[k+j],y=A[k+j+l]*g%P;
                    A[k+j]=(x+y)%P; A[k+j+l]=(x-y+P)%P;
                }
            }
        }
        if(k==-1) {ll tinv=Inv(lim); rap(i,0,lim-1) A[i]=A[i]*tinv%P;}
    }
    void MUL(ll *A,int n,ll *B,int m,ll *C){
        init(n+m); NTT(A,1); NTT(B,1); rap(i,0,lim-1) C[i]=A[i]*B[i]%P; NTT(A,-1); NTT(B,-1); NTT(C,-1);
    }
    void MUL(ll *A,int n,ll *B,int m){
        init(n+m); NTT(A,1); NTT(B,1); rap(i,0,lim-1) B[i]=B[i]*A[i]%P; NTT(A,-1); NTT(B,-1);
    }
    void INV(ll *A,ll *B,int n){
        if(n==1){B[0]=Inv(A[0]); return;} INV(A,B,(n+1)>>1);
        init(2*n); m(c,0); rap(i,0,n-1) c[i]=A[i]; NTT(c,1); NTT(B,1);
        rap(i,0,lim-1) B[i]=(2-B[i]*c[i]+P)%P*B[i]%P;
        NTT(B,-1); rap(i,n,lim-1) B[i]=0; return;
    }
    void MUL(ll *A,int *a,int n){
        int len=1; while(len<n) len<<=1; ll *f[N];
        rap(i,1,len){f[i]=new ll[lb(i)<<2]; f[i][0]=1; f[i][1]=a[i];}
        for(int l=1;l<len;l<<=1){init(2*l); for(int j=l*2;j<=len;j+=l*2){
            NTT(f[j-l],1); NTT(f[j],1);
            rap(i,0,lim-1) f[j][i]=f[j][i]*f[j-l][i]%P; NTT(f[j],-1);
        }}
        rap(i,0,n) A[i]=f[len][i];
        return;
    }
    void DEAL(ll *A,int *a,int n,int k){
        m(f,0); MUL(f,a,n); m(f1,0); rap(i,0,n-1) f1[i]=f[i+1]*(i+1)%P;
        m(fi,0); INV(f,fi,k+1); MUL(f1,n-1,fi,k,A); drap(i,k,1) A[i]=A[i-1]*inv[i]%P;
        rap(i,k+1,lim-1) A[i]=0;
        rap(i,1,k){A[i]=A[i]*i%P; if((i&1)==0) A[i]=P-A[i];}
    }
}
int n,m,K,a[N],b[N];
ll A[N],B[N],ans[N];
int main(){
    rd(n); rd(m); rap(i,1,n) rd(a[i]); rap(i,1,m) rd(b[i]); rd(K);
    Finit(max(max(n,m),K)+100);
    Poly::DEAL(A,a,n,K); Poly::DEAL(B,b,m,K);
    rap(i,1,K) A[i]=A[i]*finv[i]%P,B[i]=B[i]*finv[i]%P;
    A[0]=n; B[0]=m; Poly::MUL(A,K,B,K,ans);
    ll tinv=Inv(1ll*n*m%P); rap(i,1,K) ans[i]=ans[i]*fac[i]%P*tinv%P;
    rap(i,1,K) printf("%lld\n",ans[i]);
    return 0;
}