多项式全家桶

这是一篇用来补上一些咕掉的多项式的博客

闲着没事干，干脆从最基础的开始填好了

1.1 FFT

这是多项式全家桶最基础的一个东西，后面的多项式操作大多都建立在 \(FFT\) 的基础上

这个比较麻烦和基础，咕了

link

1.2 麦克劳林级数

在生成函数中常用到，虽然可以直接背结论

对于一个函数 \(f(x)\)，如果存在

\begin{aligned}
f(x)=F(x)=\sum a_ix^i
\end{aligned}

则称 \(F(x)\) 为 \(f(x)\) 的幂级数。

而麦克劳林级数则是一种求得 \(a_i\) 的方法

\[\begin{aligned} f(0)&=a_0\\ f'(0)&=1a_1\\ f''(0)&=1\times2 a_2\\ \vdots \\ f^{(i)}(x)&=i!a_i\\ \end{aligned} \]

那么就有:

\[\begin{aligned} f(x)&=\sum \frac{f^{(i)}(0)}{i!}x^i \end{aligned} \]

由此，可以得到一些比较特殊的式子。

举两个 \(EGF\) 中比较常用到的例子:

  1. \(ln(1-x)\)

\[\begin{aligned} &ln(1-x)'=-(1-x)^{-1}\\ &ln(1-x)''=(-1)^{3} \times (1)(1-x)^{-2}\\ \vdots\\ &ln(1-x)^{(i)}=- (i-1)!(1-x)^{-i}\\ \therefore&\ln(1-x)=-\sum^{\infty}_{i=1}\frac{x^i}{i}\\ \end{aligned} \]

  1. \(e^x\)

\[\begin{aligned} \because& \ (e^x)^{(i)}=e^x\\ \therefore& \ f(x)=e^x=\sum^{\infty}_{i=0} \frac{x^i}{i!}\\ \end{aligned} \]

1.3 多项式牛顿迭代

这是一个比较好用的多项式推到过程中的工具，它可以求出使多项式 \(G(F(x))\equiv 0 \mod x^n\) ，的零点\(F(x)\)。

大部分人学多项式模板的时候都是把这个放在后面再学，而这边则先提前讲了，因为有了牛迭之后算 \(exp\) 和 \(ln\) 都会比较方便。

考虑到当前如果已经得到了 $

\begin{aligned}
G(F_0(x))\equiv 0 \mod x^{\lceil\frac{n}{2}\rceil}
\end{aligned}
$ 时的解。

那么将 \(G(F(x))\) 在 \(F_0(x)\) 处泰勒展开，可以得到以下的式子：
\( \begin{aligned} &G(F(x))=\sum^{}_{i=0}\frac{G^{(i)(F_0(x))}}{i!}(F(x)-F_0(x))^i\\ \end{aligned} \)
对它做一些美妙推导，可以得到：

\[\begin{aligned} &\because F(x)-F_0(x) \equiv 0 \mod x^{\lceil\frac{n}{2}\rceil}\\ &\therefore G(F(x))\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x)) \mod x^{n}\\ &F(x)\equiv F_0(x)-\frac{G(F(x))}{G'(F_0(x))} \mod x^{n}\\ \end{aligned} \]

至此，可以作到倍增的方式求得最后的 \(F(x) \mod x^{n}\)。

而在大多数情况下，求解一次的复杂度是 \(O(n\log n)\)。

即 \(T(n)=T(\frac{n}{2})+O(n\log n)\)，可以得到总复杂度也是 \(O(n\log n)\) 虽然牛迭的常数经常被两只log的分治FFT吊打

2.1 多项式求逆

对于多项式h(x),求多项式使得\(f(x)*h(x) \equiv 1 \mod x^{n}\)

这是上文牛迭的一个比较基础的运用。

设 \(G(f(x))=\frac{1}{f(x)}-h(x)\)，不难发现其零点出恰好是 \(\frac{1}{f(x)}=h(x)\) 。

先特判 \(f[x^0] =\frac{1}{h[x^0]}\)。

设 \(f_0(x)\) 表示当前已经求解到的多项式，有

\[\begin{aligned} f(x)& \equiv f_0(x)-\frac{\frac{1}{f_0(x)}-h(x)}{-\frac{1}{f_0^2(x)}} \mod x^n \\ & \equiv f_0(x)(2-f_0(x)*h(x)) \mod x^n \end{aligned} \]

注意：此处的 \(h(x)\) 应当是在模 \(x^n\) 意义下的，而非 \(x^{\lceil\frac{n}{2}\rceil}\)。

2.2 多项式 \(Ln\)

求 \(f(x)\) 使得给定多项式 \(\ln h(x)=f(x)\)

这里面的一个思想也是常用的，即对于一个多项式，可以先微分再积分。

推到过程如下：

\[\begin{aligned} f(x)&=\ln h(x)\\ &=\int (\ln h(x))'\\ &=\int \frac{h'(x)}{h(x)}\\ \end{aligned} \]

多项式的求导和积分都可以 \(O(n)\) 实现，总复杂度为 \(O(n \log n)\)

2.3 多项式exp

求 \(f(x)\) 使得给定多项式 \(e^{h(x)}=f(x)\)

大致思路还是同 2.1，因此此处不再多解释，直接上推导过程。

设 \(G(f(x))= \ln f(x)-h(x)\)。

\[\begin{aligned} f(x)&\equiv f_0(x)-\frac{\ln f_0(x)-h(x)}{\frac{1}{f_0(x)}}\\ &\equiv f_0(x)(1-\ln f_0(x)+h(x)) \end{aligned} \]

总复杂度依旧为 \(O(n \log n)\)，但是由于求 \(\ln\) 的过程需要求逆，从而导致两次倍增，常数巨大 被两只log的分治FFT吊打

2.4 多项式除法

求 \(Q(x),R(x)\) 使得其满足给定多项式 \(F(x)=G(x)Q(x)+R(x)\) ,其中 \(G(x)\) 为 \(m\) 次, \(Q(x)\) 为 \(n-m\) 次, \(R(x)\) 为 \(m-1\) 次

这个需要用到一些技巧。

首先，显然有 \(F(\frac{1}{x})=G(\frac{1}{x})Q(\frac{1}{x})+R(\frac{1}{x})\) 成立。

对于 \(F(\frac{1}{x})\) ，如果我们让它乘上一个 \(x^n\) ，会发现其恰好是 \(F(x)\) 翻转之后的结果，即 \(x^{n}F(\frac{1}{x})=revF(x)\) 。

因此，对上式左右两边都乘上一个 \(x^{n}\) ，可以得到：

\[\begin{aligned} x^{n}F(x)&=x^{n}revQ(x)revG(x)+x^{n}revR(x)\\ &=x^{n-m}revQ(x)x^{m}revG(x)+x^{n-m+1}x^{m-1}revR(x)\\ revF(x)&=revQ(x)revG(x)+x^{n-m+1}revR(x)\\ \end{aligned} \]

注意到 \(Q(x)\) 的最高项是 \(n-m\)，因此 \(revQ(x)\) 的最高项也是 \(n-m\) ，其在模 \(x^{n-m+1}\) 的情况下仍是它本身。

那么可以让这个式子模上一个 \(x^{n-m+1}\)，来去掉 \(revR(x)\) 的影响。

\[\begin{aligned} revF(x)&=revQ(x)revG(x) \mod x^{n-m+1}\\ revQ(x)&=revF(x)revG'(x) \mod x^{n-m+1} \end{aligned} \]

这样就可以求出 \(revQ(x)\) ，也就是求出了 \(Q(x)\) 。
然后根据 \(R(x)=F(x)-G(x)Q(x)\) 也可以求出 \(R(x)\)。

2.5 多项式快速幂

这个一般有两种做法。

1.暴力快速幂

可以仿照普通快速幂的样子，只不过把乘法改为多项式乘法即可。

时间复杂度 \(O(n\log n \log k)\) ，效率比较低，但不失为一种好的思路。比如在线性递推中，就需要用到这个思想。

不清楚对于多项式能否使用费马小定理来解决 \(k\) 比较大的时候的情况。

2. \(\exp+ \ln\)

根据高中数学书上的内容， \(\ln x^{k}=k \ln x\)。

这个在 \(x\) 是一个多项式的时候同样适用。

因此，可以得到 \(F(x)^k=\exp(k \ln F(x))\)。

时间复杂度 \(O(n \log n)\)。

而且当 \(k\) 是个比较大的数（需要取模）的时候，这个方法也同样适用。

不过不是让 \(k\) 对 \(P-1\) 取模，而是对 \(P\) 取模。

给出一个不太严谨的理解方式。

\(e^{kx}=\sum^{\infty}_{i=0}\frac{k^i x^i}{i!}\),因此 \(k\) 模 \(P\) 是对的。

2.6 板子

到这里，已经差不多讲完了多项式的几个常用的板子，下面放出本人的模板不开O2经常会被卡。

数组版 
namespace Poly{
    const ll mod=998244353;
    ll ksm(ll d,ll t){
        ll res=1;
        for(;t;t>>=1){
            if(t&1) res=res*d%mod;
            d=d*d%mod;
        }
        return res;
    }
    il ll dec(ll x,ll y) {return x>=y?x-y:x-y+mod;}
    il ll add(ll x,ll y) {return x+y<mod?x+y:x+y-mod;}
    const int N=262144;
    vector<ll> w[23];
    ll fac[N],ifac[N],inv[N];
    il ll C(ll x,ll y){
        if(x<y) return 0;
        return fac[x]*ifac[y]%mod*ifac[x-y]%mod;
    }
    void init(){
        int n=N-1;
        fac[0]=1;
        for(ri i=1;i<=n;++i) fac[i]=i*fac[i-1]%mod;
        ifac[n]=ksm(fac[n],mod-2);
        for(ri i=n-1;~i;--i) ifac[i]=ifac[i+1]*(i+1)%mod;
        for(ri i=n;i;--i) inv[i]=ifac[i]*fac[i-1]%mod;
        inv[0]=1;
        int d=log(N)/log(2)+0.5;
        for(ri i=1;i<=d;++i){
            w[i].resize(1<<i);
            w[i][0]=1,w[i][1]=ksm(3,(mod-1)>>i);
            for(ri j=2;j<(1<<i);++j) w[i][j]=w[i][j-1]*w[i][1]%mod;
        }
    }
    int r[N];
    void DFT(int limit,ll *a,int flag){
        for(ri i=1;i<limit;++i) if(r[i]<i) swap(a[i],a[r[i]]);
        for(ri l=1,t=1;l<limit;l<<=1,++t){
            for(ri i=0;i<limit;i+=l<<1){
                ll *W=&w[t][0];
                for(ri j=0;j<l;++j){
                    ll tmp=a[i+j+l]*(*W++)%mod;
                    a[i+j+l]=dec(a[i+j],tmp);
                    a[i+j]=add(a[i+j],tmp);
                }
            }
        }
        if(flag==-1){
            reverse(a+1,a+limit);
            ll inv=ksm(limit,mod-2);
            for(ri i=0;i<limit;++i) a[i]=a[i]*inv%mod;
        }
    }
    ll _F[N];
    void rev(int limit,int ws){
        for(ri i=0;i<limit;i+=2){
            r[i]=r[i>>1]>>1;
            r[i|1]=(r[i>>1]>>1)|(1<<ws-1);
        }
    }
    void NTT(int n,ll *f,int m,ll *g,int lim=0){
        if(!lim) lim=n+m;
        int limit=1,ws=0;
        for(;limit<=n+m;++ws,limit<<=1);
        rev(limit,ws);
        for(ri i=0;i<=m;++i) _F[i]=g[i];
        for(ri i=m+1;i<limit;++i) _F[i]=0;
        for(ri i=n+1;i<limit;++i) f[i]=0;
        DFT(limit,f,1),DFT(limit,_F,1);
        for(ri i=0;i<limit;++i) f[i]=f[i]*_F[i]%mod;
        DFT(limit,f,-1);
        for(ri i=lim+1;i<limit;++i) f[i]=0;
    }
    void Inv(int n,ll *h,ll *f){
        memset(_F,0,sizeof(_F));
        f[0]=ksm(h[0],mod-2);
        for(ri t=1,l=2;1;l<<=1,++t){
            for(ri i=0;i<l;++i) _F[i]=h[i];
            rev(l<<1,t+1);
            DFT(l<<1,f,1),DFT(l<<1,_F,1);
            for(ri i=0;i<(l<<1);++i) f[i]=(2*f[i]-f[i]*f[i]%mod*_F[i]%mod+mod)%mod;
            DFT(l<<1,f,-1);
            if(l>n){
                for(ri i=n+1;i<(l<<1);++i) f[i]=0;
                break;
            }
            for(ri i=l;i<(l<<1);++i) f[i]=0;
        }
    }
    il void der(int n,ll *f){
        for(ri i=0;i<n;++i) 
            f[i]=f[i+1]*(i+1)%mod;
        f[n]=0;
    }
    il void Int(int n,ll *f){
        for(ri i=n;i;--i) 
            f[i]=f[i-1]*inv[i]%mod;
        f[0]=0;
    }
    ll _G[N];
    il void Ln(int n,ll *h,ll *f){
        memset(_G,0,sizeof(_G));
        for(ri i=0;i<=n;++i) f[i]=h[i];
        Inv(n,h,_G);
        der(n,f);
        NTT(n,f,n,_G);
        Int(n,f);
    }
    ll exp_f[N];
    il void exp(int n,ll *h,ll *f){
        f[0]=1;
        for(ri l=2,t=1;1;l<<=1,t++){
            memset(exp_f,0,sizeof(exp_f));
            Ln(l-1,f,exp_f);
            for(ri i=0;i<l;++i) exp_f[i]=(-exp_f[i]+h[i]+mod)%mod;
            exp_f[0]++;
            NTT(l>>1,f,l,exp_f);
            for(ri i=l;i<(l<<1);++i) f[i]=0;
            if(l>n){
                for(ri i=n+1;i<l;++i) f[i]=0;
                break;
            }
        }
    }
    ll ksm_f[N];
    il void polyksm(int n,ll *f,ll t){
        memset(ksm_f,0,sizeof(ksm_f));
        for(ri i=0;i<=n;++i) ksm_f[i]=f[i];
        Ln(n,ksm_f,f);
        for(ri i=0;i<=n;++i) ksm_f[i]=f[i]*t%mod,f[i]=0;
        exp(n,ksm_f,f);
    }
    il void div(int n,ll *f,int m,ll *g,ll *Q,ll *R){
        reverse(f,f+n+1),reverse(g,g+m+1);
        Inv(n-m,g,Q),NTT(n-m,Q,n-m,f,n-m);
        reverse(f,f+n+1),reverse(g,g+m+1),reverse(Q,Q+n-m+1);
        for(ri i=0;i<=m;++i) R[i]=g[i];NTT(m,R,n-m,Q);
        for(ri i=0;i<=n;++i) R[i]=dec(f[i],R[i]);
    }
}
没有写快速幂和除法的更慢的vector版本，不过用起来会方便一些 
const ll mod=998244353;
namespace Poly{
    il ll ksm(ll d,ll t){
        ll res=1;
        for(;t;t>>=1,d=d*d%mod)
            if(t&1) res=res*d%mod;
        return res;
    }
    const int N=2097152;
    vector<ll> w[23];
    ll fac[N],ifac[N],inv[N];
    il ll C(ll x,ll y){
        if(x<y) return 0;
        return fac[x]*ifac[y]%mod*ifac[x-y]%mod;
    }
    void init(){
        int n=N-1;
        fac[0]=1;
        for(ri i=1;i<=n;++i) fac[i]=i*fac[i-1]%mod;
        ifac[n]=ksm(fac[n],mod-2);
        for(ri i=n-1;~i;--i) ifac[i]=ifac[i+1]*(i+1)%mod;
        inv[0]=1;
        for(ri i=1;i<=n;++i) inv[i]=fac[i-1]*ifac[i]%mod;
        int d=log(N)/log(2)+0.5;
        for(ri t=1;t<=d;++t){
            w[t].resize(1<<t);
            w[t][0]=1,w[t][1]=ksm(3,(mod-1)>>t);
            for(ri i=2;i<(1<<t);++i) w[t][i]=w[t][i-1]*w[t][1]%mod;
        }
    }
    int r[N];
    il void rev(int limit,int ws){
        for(ri i=0;i<limit;i+=2){
            r[i]=r[i>>1]>>1;
            r[i|1]=(r[i>>1]>>1)|(1<<ws-1);
        }
    }
    il void DFT(int limit,vector<ll> &f,int flag){
        for(ri i=1;i<limit;++i) 
            if(i>r[i]) swap(f[i],f[r[i]]);
        for(ri t=1,l=1;l<limit;++t,l<<=1){
            for(ri i=0;i<limit;i+=l<<1){
                ll *W=&w[t][0];
                for(ri j=0;j<l;++j){
                    ll tmp=f[i+j+l]*(*W++)%mod;
                    f[i+j+l]=f[i+j]-tmp;
                    f[i+j]=f[i+j]+tmp;
                    if(f[i+j+l]<0) f[i+j+l]+=mod;
                    if(f[i+j]>=mod) f[i+j]-=mod;
                }
            }
        }
        if(flag==-1){
            ll inv=ksm(limit,mod-2);
            reverse(f.begin()+1,f.begin()+limit);
            for(ri i=0;i<limit;++i) f[i]=f[i]*inv%mod;
        }
    }
    il vector<ll> NTT(int n,vector<ll> f,int m,vector<ll> g){
        int limit=1,ws=0;
        for(;limit<=n+m;limit<<=1,++ws);
        f.resize(limit),g.resize(limit);
        f.erase(f.begin()+limit,f.end());
        g.erase(g.begin()+limit,g.end());
        rev(limit,ws);
        DFT(limit,f,1),DFT(limit,g,1);
        for(ri i=0;i<limit;++i) f[i]=f[i]*g[i]%mod;
        DFT(limit,f,-1);
        f.erase(f.begin()+n+m+1,f.end());
        return f;
    }
    il vector<ll> Inv(int n,vector<ll> &h){
        vector<ll> f,g;
        f.resize(n<<2),g.resize(n<<2);
        f[0]=ksm(h[0],mod-2);
        for(ri t=1,l=2;1;l<<=1,++t){
            for(ri i=0;i<l;++i) g[i]=h[i];
            rev(l<<1,t+1);
            DFT(l<<1,f,1),DFT(l<<1,g,1);
            for(ri i=0;i<(l<<1);++i) f[i]=(2*f[i]-f[i]*f[i]%mod*g[i]%mod+mod)%mod;
            DFT(l<<1,f,-1);
            for(ri i=l;i<(l<<1);++i) f[i]=0;
            if(l>n) break;
        }
        f.erase(f.begin()+n+1,f.end());
        return f;                
    }
    il vector<ll> der(int n,vector<ll> f){
        for(ri i=0;i<n;++i) 
            f[i]=f[i+1]*(i+1)%mod;
        f[n]=0;
        return f;
    }
    il vector<ll> Int(int n,vector<ll> f){
        for(ri i=n;i;--i) f[i]=f[i-1]*inv[i]%mod;
        f[0]=0;
        return f;
    }
    il vector<ll> Ln(int n,vector<ll> f){
        vector<ll> g=Inv(n,f);
        f=der(n,f);
        f=NTT(n,f,n,g);
        f.erase(f.begin()+n+1,f.end());
        f=Int(n,f);
        return f;
    }
    il vector<ll> exp(int n,vector<ll> &h){
        vector<ll> f,g,F;
        int limit=1;
        f.resize(n<<2);
        f[0]=1;
        for(ri l=2,t=1;1;l<<=1,t++){
            g.resize(l);
            for(ri i=0;i<l;++i) g[i]=f[i];
            g=Ln(l-1,g);
            for(ri i=0;i<l;++i) g[i]=(-g[i]+h[i]+mod)%mod;
            g[0]++;
            F=NTT(l>>1,f,l,g);
            for(ri i=l>>1;i<l;++i) f[i]=F[i];
            if(l>n) break;
        }
        f.erase(f.begin()+n+1,f.end());
        return f;
    }
}

3.1 多项式多点求值

给你一个 \(n\) 次多项式 \(f(x)\) 和一个长为 \(m\) 的序列 \(a\) 让你求出 \(x\) 分别为 \(a_{i}\) 时的值。

方法1：我会暴力！

直接秦九昭展开。

时间复杂度 \(O(nm)\) ，空间复杂度 \(O(n+m)\)。

方法2：我会多项式除法！

回忆 2.4 中讲的多项式除法，还记得这样一个式子吗？

\(\begin{aligned}F(x)&=Q(x)G(x)+R(x)\end{aligned}\)

如果构造 \(G(x)=x-a_{i}\) ，不是就有 \(F(a_{i})=Q(a_{i})\times 0 +R(a_{i})=R(a_{i})\) 了？

于是，就得到了一个 O(nm log m) 的优秀做法

首先，取模有一个非常显然的性质： \(x \mod a=(x \mod ab)\mod a,b \geq 1\)

这个在多项式意义下也是成立的，因此考虑分治。

具体地说，假设当前需要得到 \(F(x) \mod (x-a_{i}),i\in[l,r]\)，可以先求出\(F(x) \mod \prod^{r}_{i=l}(x-a_{i})\)。

然后可以把得到的结果继续代入到区间 \([l,mid],[mid+1,r]\) 中去做。

因为取模的性质，每次取模之后的 $F(x) $ 项数都会是 \(r-l+1\)。

那么有 \(T(n)=2T(\frac{n}{2})+O(n \log n)\) ，复杂度是 \(O(n \log^2 n)\) 的。

而求 \(\prod^{r}_{i=l} (x-a_{i})\) 的过程也可以使用分治，复杂度和上面的分析方法一样，做的过程中可以拿vector或者指针记录下区间 \([l,r]\) 的结果。

因此这样，就可以得到一个时间复杂度 \(O(n \log^2 n)\)，空间复杂度 \(O(n \log n)\) 的较为优秀的做法。

然后就会发现跑一组n,m=64000的数据要1m

多项式取模的常数非常大，而当区间比较小的时候会频繁调用取模。

因此，可以当区间长度小于某个阈值的时候使用常数很小的暴力取模，可以起到非常大的优化。

Code 
#include<cstdio>
#include<vector>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
#define il inline
#define ri register int
#define ll long long
#define ui unsigned int
il ll read(){
    bool f=true;ll x=0;
    register char ch=getchar();
    while(ch<'0'||ch>'9') {if(ch=='-') f=false;ch=getchar();}
    while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
    if(f) return x;
    return ~(--x);
}
il void write(const ll &x){if(x>9) write(x/10);putchar(x%10+'0');}
il void print(const ll &x) {x<0?putchar('-'),write(~(x-1)):write(x);putchar('\n');}
il ll max(const ll &a,const ll &b){return a>b?a:b;}
il ll min(const ll &a,const ll &b){return a<b?a:b;}
namespace Poly{
    const ll mod=998244353;
    ll ksm(ll d,ll t){
        ll res=1;
        for(;t;t>>=1){
            if(t&1) res=res*d%mod;
            d=d*d%mod;
        }
        return res;
    }
    il ll dec(ll x,ll y) {return x>=y?x-y:x-y+mod;}
    il ll add(ll x,ll y) {return x+y<mod?x+y:x+y-mod;}
    const int N=262144;
    vector<ll> w[23];
    ll fac[N],ifac[N],inv[N];
    il ll C(ll x,ll y){
        if(x<y) return 0;
        return fac[x]*ifac[y]%mod*ifac[x-y]%mod;
    }
    void init(){
        int n=N-1;
        fac[0]=1;
        for(ri i=1;i<=n;++i) fac[i]=i*fac[i-1]%mod;
        ifac[n]=ksm(fac[n],mod-2);
        for(ri i=n-1;~i;--i) ifac[i]=ifac[i+1]*(i+1)%mod;
        for(ri i=n;i;--i) inv[i]=ifac[i]*fac[i-1]%mod;
        inv[0]=1;
        int d=log(N)/log(2)+0.5;
        for(ri i=1;i<=d;++i){
            w[i].resize(1<<i);
            w[i][0]=1,w[i][1]=ksm(3,(mod-1)>>i);
            for(ri j=2;j<(1<<i);++j) w[i][j]=w[i][j-1]*w[i][1]%mod;
        }
    }
    int r[N];
    void DFT(int limit,ll *a,int flag){
        for(ri i=1;i<limit;++i) if(r[i]<i) swap(a[i],a[r[i]]);
        for(ri l=1,t=1;l<limit;l<<=1,++t){
            for(ri i=0;i<limit;i+=l<<1){
                ll *W=&w[t][0];
                for(ri j=0;j<l;++j){
                    ll tmp=a[i+j+l]*(*W++)%mod;
                    a[i+j+l]=dec(a[i+j],tmp);
                    a[i+j]=add(a[i+j],tmp);
                }
            }
        }
        if(flag==-1){
            reverse(a+1,a+limit);
            ll inv=ksm(limit,mod-2);
            for(ri i=0;i<limit;++i) a[i]=a[i]*inv%mod;
        }
    }
    ll _F[N];
    void rev(int limit,int ws){
        for(ri i=0;i<limit;i+=2){
            r[i]=r[i>>1]>>1;
            r[i|1]=(r[i>>1]>>1)|(1<<ws-1);
        }
    }
    void NTT(int n,ll *f,int m,ll *g,int lim=0){
        if(!lim) lim=n+m;
        int limit=1,ws=0;
        for(;limit<=n+m;++ws,limit<<=1);
        rev(limit,ws);
        for(ri i=0;i<=m;++i) _F[i]=g[i];
        for(ri i=m+1;i<limit;++i) _F[i]=0;
        for(ri i=n+1;i<limit;++i) f[i]=0;
        DFT(limit,f,1),DFT(limit,_F,1);
        for(ri i=0;i<limit;++i) f[i]=f[i]*_F[i]%mod;
        DFT(limit,f,-1);
        for(ri i=lim+1;i<limit;++i) f[i]=0;
    }
    void Inv(int n,ll *h,ll *f){
        memset(_F,0,sizeof(_F));
        f[0]=ksm(h[0],mod-2);
        for(ri t=1,l=2;1;l<<=1,++t){
            for(ri i=0;i<l;++i) _F[i]=h[i];
            rev(l<<1,t+1);
            DFT(l<<1,f,1),DFT(l<<1,_F,1);
            for(ri i=0;i<(l<<1);++i) f[i]=(2*f[i]-f[i]*f[i]%mod*_F[i]%mod+mod)%mod;
            DFT(l<<1,f,-1);
            if(l>n){
                for(ri i=n+1;i<(l<<1);++i) f[i]=0;
                break;
            }
            for(ri i=l;i<(l<<1);++i) f[i]=0;
        }
    }
    il void div(int n,ll *f,int m,ll *g,ll *Q,ll *R){
        reverse(f,f+n+1),reverse(g,g+m+1);
        Inv(n-m,g,Q),NTT(n-m,Q,n-m,f,n-m);
        reverse(f,f+n+1),reverse(g,g+m+1),reverse(Q,Q+n-m+1);
        for(ri i=0;i<=m;++i) R[i]=g[i];NTT(m,R,n-m,Q);
        for(ri i=0;i<=n;++i) R[i]=dec(f[i],R[i]);
    }
    ll Mod_Q[N],Mod_f[N];
    void Mod_short(int n,ll *f,int m,ll *g){
        for(ri i=n;i>=m;--i){
            if(f[i]){
                ll x=f[i]*ksm(g[m],mod-2)%mod;
                for(ri j=0;j<=m;++j){
                    f[i-j]=dec(f[i-j],g[m-j]*x%mod);
                }
            }
        }
    }
    il void Mod(int n,ll *f,int m,ll *g){
        if(n-m<=400){
            Mod_short(n,f,m,g);
            return;
        }
        // if(n<m) return;
        for(ri i=0;i<=n;++i) Mod_Q[i]=0;
        for(ri i=0;i<=n;++i) Mod_f[i]=f[i],f[i]=0;
        reverse(Mod_f,Mod_f+n+1),reverse(g,g+m+1);
        Inv(n-m,g,Mod_Q),NTT(n-m,Mod_Q,n-m,Mod_f,n-m);
        reverse(Mod_f,Mod_f+n+1),reverse(g,g+m+1),reverse(Mod_Q,Mod_Q+n-m+1);
        for(ri i=0;i<=m;++i) f[i]=g[i];NTT(m,f,n-m,Mod_Q);
        for(ri i=0;i<=n;++i) f[i]=dec(Mod_f[i],f[i]);
    }
}
using namespace Poly;
/*
两次分治，第一次求乘起来的，第二次求取模结果
*/
ll G[N*32],*g[N<<4],a[N],res[N];
int lst;
void apply(int id,int len){//动态申请空间
    g[id]=G+lst;
    lst+=len;
}
int cl[N],cr[N],cnt;
int solve_1(int l,int r){
    int now=++cnt;
    cl[cnt]=l,cr[cnt]=r;
    apply(cnt,r-l+1+1);
    if(l==r){
        g[cnt][0]=mod-a[l];
        g[cnt][1]=1;
        return now;
    }
    int mid=(l+r)>>1;
    int L=solve_1(l,mid),R=solve_1(mid+1,r);
    for(ri i=0;i<=mid-l+1;++i) res[i]=g[L][i];
    NTT(mid-l+1,res,r-mid,g[R]);
    for(ri i=0;i<=r-l+1;++i) g[now][i]=res[i];
    return now;
}
ll ans[N];
void solve_2(int l,int r,vector<ll> &f){
    int now=++cnt;
    if(r-l+1<f.size()){
        for(ri i=0;i<f.size();++i) res[i]=f[i];
        res[f.size()]=0;
        Mod(f.size(),res,r-l+1,g[now]);
        for(ri i=0;i<r-l+1;++i) f[i]=res[i];
        f.erase(f.begin()+r-l+1,f.end());
    }
    if(l==r){
        ans[l]=f[0];
        return;
    }
    int mid=(l+r)>>1;
    vector<ll> rf=f;
    solve_2(l,mid,f),solve_2(mid+1,r,rf);
}
vector<ll> f;
int n,m;
int main(){
    // freopen("P5050_6.in","r",stdin);
    // freopen("1.out","w",stdout);
    init();
    n=read(),m=read();
    f.resize(n+1);
    for(ri i=0;i<=n;++i) f[i]=read();
    for(ri i=1;i<=m;++i) a[i]=read();
    m=max(n,m);
    solve_1(1,m),cnt=0;
    solve_2(1,m,f);
    for(ri i=1;i<=m;++i) print(ans[i]);
    return 0;
}

3.2 普通幂转下降幂

给你一个普通幂多项式 \(F(x)\) ，求出其下降幂多项式 \(G(x)\) 使得 \(G(x)=F(x)\) 。

首先，要解释什么是下降幂。

对于 \(x^{\underline{n}}\) ，满足递推式 \(x^{\underline{n}}= \left \{ \begin{aligned} &1,n=0\\ &x\times (x-1)^{\underline{n-1}},n>0\\ \end{aligned} \right . \)

换种写法就是 \(x^{\underline{n}}=\frac{x!}{(x-n)!}\) 。

回忆推导第二类斯特林数的时候，那个经典的反演式子。

\[\begin{aligned} m^n&=\sum^{n}_{i=0} \left \{ \matrix{n\\i}\right \}C^{m}_{i}i!\\ &=\sum^{n}_{i=0} \left \{ \matrix{n\\i}\right \}\frac{m!}{(m-i)!}\\ &=\sum^{n}_{i=0} \left \{ \matrix{n\\i}\right \}m^{\underline{i}}\\ \end{aligned} \]

尝试把 \(m\) 换成 \(x\)，那么有：

\[\begin{aligned} x^{n}&=\sum_{i=0}^{n}\left \{ \matrix{n\\i}\right \}x^{\underline{i}}\\ F(x)&=\sum^{n}_{i=0}f_{i}\sum^{i}_{j=0}\left \{ \matrix{i\\j}\right \}x^{\underline{j}}\\ &=\sum^{n}_{i=0}x^{\underline{i}}\sum^{n}_{j=i}\left \{ \matrix{j\\i}\right \}f_{j}\\ \therefore a_{i}&=\sum^{n}_{j=i}\left \{ \matrix{j\\i}\right \}f_{j}\\ &=\sum^{n}_{j=0}\left \{ \matrix{j\\i}\right \}f_{j}\\ &=\sum^{n}_{j=0}f_{j}\sum^{i}_{k=0}\frac{(-1)^{i-k}k^j}{k!(i-k)!}\\ &=\sum^{i}_{k=0}\frac{(-1)^{i-k}}{k!(i-k)!}\sum^{n}_{j=0}f_{j}k^j\\ &=\sum^{i}_{k=0}\frac{(-1)^{i-k}}{k!(i-k)!}g_{k}\\ \end{aligned} \]

右边部分显然是一个卷积的形式，那么此时瓶颈在于如何得到 \(g_{k}\)

这里是多项式多点求值啊......

那么在经过多项式多点求值之后，把\(g(x)\) 和左边那一陀卷起来，就得到答案啦！

非常离谱的是这些东西是我在尝试推下降幂多项式卷积的时候得到的东西

3.3 下降幂多项式卷积

给两个下降幂多项式 \(F(x)\) 与 \(G(x)\) ，求它们的卷积。

按照之前 \(FFT\) 的思路，只要求出了 \(F(x)\) 和 \(G(x)\) 的点值表达式，把它们乘起来再转回系数表达式。

设\(f(x)\) 是 \(F(x)\)的点值的生成函数。

\[\begin{aligned} f[m]&=\sum^{m}_{i=0}\frac{m!}{(m-i)!}a_{i}\\ \frac{f[m]}{m!}&=\sum^{m}_{i=0}\frac{a_{i}}{(m-i)!} \end{aligned} \]

发现改为 \(EGF\) 形式之后就是 \(F(x)e^{x}\)

那么从点值变回来也只要卷上一个 \(e^{-x}\) 就好了。

Code 
#include<bits/stdc++.h>
using namespace std;
#define il inline
#define ri register int
#define ll long long
#define ui unsigned int
il ll read(){
    bool f=true;ll x=0;
    register char ch=getchar();
    while(ch<'0'||ch>'9') {if(ch=='-') f=false;ch=getchar();}
    while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
    if(f) return x;
    return ~(--x);
}
il void write(const ll &x){if(x>9) write(x/10);putchar(x%10+'0');}
il void print(const ll &x) {x<0?putchar('-'),write(~(x-1)):write(x);putchar(' ');}
il ll max(const ll &a,const ll &b){return a>b?a:b;}
il ll min(const ll &a,const ll &b){return a<b?a:b;}
namespace Poly{
    const ll mod=998244353;
    ll ksm(ll d,ll t){
        ll res=1;
        for(;t;t>>=1){
            if(t&1) res=res*d%mod;
            d=d*d%mod;
        }
        return res;
    }
    il ll dec(ll x,ll y) {return x>=y?x-y:x-y+mod;}
    il ll add(ll x,ll y) {return x+y<mod?x+y:x+y-mod;}
    const int N=262144<<1;
    vector<ll> w[23];
    ll fac[N],ifac[N],inv[N];
    il ll C(ll x,ll y){
        if(x<y) return 0;
        return fac[x]*ifac[y]%mod*ifac[x-y]%mod;
    }
    void init(){
        int n=N-1;
        fac[0]=1;
        for(ri i=1;i<=n;++i) fac[i]=i*fac[i-1]%mod;
        ifac[n]=ksm(fac[n],mod-2);
        for(ri i=n-1;~i;--i) ifac[i]=ifac[i+1]*(i+1)%mod;
        for(ri i=n;i;--i) inv[i]=ifac[i]*fac[i-1]%mod;
        inv[0]=1;
        int d=log(N)/log(2)+0.5;
        for(ri i=1;i<=d;++i){
            w[i].resize(1<<i);
            w[i][0]=1,w[i][1]=ksm(3,(mod-1)>>i);
            for(ri j=2;j<(1<<i);++j) w[i][j]=w[i][j-1]*w[i][1]%mod;
        }
    }
    int r[N];
    void DFT(int limit,ll *a,int flag){
        for(ri i=1;i<limit;++i) if(r[i]<i) swap(a[i],a[r[i]]);
        for(ri l=1,t=1;l<limit;l<<=1,++t){
            for(ri i=0;i<limit;i+=l<<1){
                ll *W=&w[t][0];
                for(ri j=0;j<l;++j){
                    ll tmp=a[i+j+l]*(*W++)%mod;
                    a[i+j+l]=dec(a[i+j],tmp);
                    a[i+j]=add(a[i+j],tmp);
                }
            }
        }
        if(flag==-1){
            reverse(a+1,a+limit);
            ll inv=ksm(limit,mod-2);
            for(ri i=0;i<limit;++i) a[i]=a[i]*inv%mod;
        }
    }
    ll _F[N];
    void rev(int limit,int ws){
        for(ri i=0;i<limit;i+=2){
            r[i]=r[i>>1]>>1;
            r[i|1]=(r[i>>1]>>1)|(1<<ws-1);
        }
    }
}
using namespace Poly;
int n,m;
ll f[N],g[N],e[N];
int main(){
    init();
    n=read(),m=read();
    for(ri i=0;i<=n;++i) f[i]=read();
    for(ri i=0;i<=m;++i) g[i]=read();
    for(ri i=0;i<=n+m;++i) e[i]=ifac[i];
    int limit=1,ws=0;
    for(;limit<=(n+m)*2;limit<<=1,++ws);
    rev(limit,ws);
    DFT(limit,e,1),DFT(limit,f,1),DFT(limit,g,1);
    for(ri i=0;i<limit;++i) f[i]=f[i]*e[i]%mod,g[i]=g[i]*e[i]%mod,e[i]=0;
    DFT(limit,f,-1),DFT(limit,g,-1);
    for(ri i=n+m+1;i<limit;++i) f[i]=0;
    for(ri i=0;i<=n+m;++i){
        f[i]=f[i]*g[i]%mod*fac[i]%mod;
        if(i&1) e[i]=mod-ifac[i];
        else e[i]=ifac[i];
    }
    DFT(limit,f,1),DFT(limit,e,1);
    for(ri i=0;i<limit;++i) f[i]=f[i]*e[i]%mod;
    DFT(limit,f,-1);
    for(ri i=n+m+1;i<limit;++i) f[i]=0;
    for(ri i=0;i<=n+m;++i) print(f[i]);
    return 0;
}
posted @ 2021-04-30 10:40  krimson  阅读(295)  评论(1)  编辑  收藏  举报