【多项式】advanced+生成函数

前言

发现多项式的坑挺深，一文堆不完。所以另开新坑。
基础多项式算法

分治FFT

洛谷板子

在上一篇文章中，我们用生成函数的思路以及多项式求逆的技术解决了这道题。
现在我们用CDQ分治的思路以及多项式乘法的技术来解决这道题。
首先，我们发现$f[i]$是由$f[0]$~$f[i-1]$与$g$卷积而来。考虑暴力（不用fft）计算，时间复杂度为$O(n^2)$。
把前几项$f[i]$的答案列出，观察：
$f_1=f_0{g}_1$
$f_2=f_0{g}_2+f_1{g}_1$
$f_3=f_0{g}_3+f_1{g}_2+f_2{g}_1$
考虑CDQ分治。
把当前处理区间$f_l$~$f_r$，设$mid=\lfloor \frac{l+r}{2}\rfloor$。
先递归计算$f_l$~$f_{mid}$，接着考虑前半部分对后半部分的贡献。
设对$f_x(mid+1\le x\le r)$的贡献为$w_x$，则有
$$w_x=\sum _{i=l}^{mid}{f}_i\times g_{x-i}$$
然后你发现这就是个卷积！
那就可以用$O(nlogn)$的时间内求出左边区间对于右边区间每个位置的贡献！
然后如何具体处理呢？
$$\sum _{i=l}^{mid}{f}_i\times g_{x-i}$$
设$f_i=0(mid+1\le i\le r)$，
$$\sum _{i=l}^x{f}_i\times g_{x-i}$$
$$\sum _{x=0}^{x-l}{f}_{i+l}{g}_{x-l-i}$$
设$a_i=f_{i+l},b_i=g_i$。
原式为
$$\sum _{x=0}^{x-l}{a}_i{b}_{x-l-i}$$
于是就能算了。
就是长度为该区间长度(r-l)/2的两个多项式相乘。得到一个(r-l)的多项式。但是我们只取后面的部分(x>mid)。

// 分治fft
// CDQ分治思路:即左区间每个值对右区间每个值
// 的贡献可以O(n)或此处:O(nlogn)计算出来
// 为什么可以O(nlogn)计算出来?因为借用了卷积形式

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

const int MAXN = 1e5+5;
inline int read(){int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch = getchar();}while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}return x*f;}

const int P = 998244353;

inline int qpow(int x,int y){int res(1);while(y){if(y&1) res = 1ll*res*x%P;x = 1ll*x*x%P;y>>=1;}return res;}

int rev[MAXN<<2];

void change(int y[],int len){
    for(int i=0;i<len;++i){
        rev[i] = rev[i>>1]>>1;
        if(i&1) rev[i] |= len>>1;
    }
    for(int i=0;i<len;++i)
        if(i<rev[i]) swap(y[i],y[rev[i]]);
}

void ntt(int y[],int len,int on){
    change(y,len);
    for(int h=2;h<=len;h<<=1){
        int gn = qpow(3,(P-1)/h);
        for(int j=0;j<len;j+=h){
            int g = 1;
            for(int k=j;k<j+h/2;++k){
                int u = y[k];
                int t = 1ll*g*y[k+h/2]%P;
                y[k] = (u+t)%P;
                y[k+h/2] = (u-t+P)%P;
                g = 1ll*g*gn%P;
            }
        }
    }
    if(on==-1) {
        reverse(y+1,y+len);
        register int inv = qpow(len,P-2);
        for(int i=0;i<len;++i) y[i] = 1ll*y[i]*inv%P;
    }
}

int f[MAXN<<2],g[MAXN<<2];
int a[MAXN<<2],b[MAXN<<2]; // temp

void solve(int l,int r){
    // 分治fft:
    if(l+1>=r) return;
    int len = (r-l)>>1; // 半长
    int mid = len+l;
    solve(l,mid); // 递归计算左区间
    memset(a+len,0,sizeof(int)*len);
    memcpy(a,f+l,sizeof(int)*len); // a中储存的是?
    memcpy(b,g,sizeof(int)*(r-l)); // b中储存的是?
    ntt(a,r-l,1),ntt(b,r-l,1);
    for(int i=0;i<r-l;++i) a[i]=1ll*a[i]*b[i]%P;
    ntt(a,r-l,-1);
    for(int i=len;i<r-l;++i) f[l+i]=(1ll*f[l+i]+a[i])%P;
    solve(mid,r); // 递归计算右区间
}

int main(){
    int n=read();
    for(int i=1;i<n;++i) g[i]=read();
    f[0]=1;
    int len=1;
    while(n>len) len<<=1;
    solve(0,len);
    for(int i=0;i<n;++i)
        printf("%d%c",f[i],i==n-1?'\n':' ');
}

多项式对数函数

洛谷板子
不用动脑子的推导。
推导：

// 多项式对数函数

#include<bits/stdc++.h>
using namespace std;

const int MAXN = 1e5+5;
inline int read(){int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch = getchar();}while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}return x*f;}

const int P = 998244353;

inline int qpow(int x,int y){int res(1);while(y){if(y&1) res = 1ll*res*x%P;x = 1ll*x*x%P;y>>=1;}return res;}

int rev[MAXN<<2];

void change(int y[],int len){
    for(int i=0;i<len;++i){
        rev[i] = rev[i>>1]>>1;
        if(i&1) rev[i] |= len>>1;
    }
    for(int i=0;i<len;++i)
        if(i<rev[i]) swap(y[i],y[rev[i]]);
}

void ntt(int y[],int len,int on){
    change(y,len);
    for(int h=2;h<=len;h<<=1){
        int gn = qpow(3,(P-1)/h);
        for(int j=0;j<len;j+=h){
            int g = 1;
            for(int k=j;k<j+h/2;++k){
                int u = y[k];
                int t = 1ll*g*y[k+h/2]%P;
                y[k] = (u+t)%P;
                y[k+h/2] = (u-t+P)%P;
                g = 1ll*g*gn%P;
            }
        }
    }
    if(on==-1) {
        reverse(y+1,y+len);
        register int inv = qpow(len,P-2);
        for(int i=0;i<len;++i) y[i] = 1ll*y[i]*inv%P;
    }
}

static int inv_t[MAXN<<2];

void get_inv(int f[],int g[],int n){
    // 递归:把两个n/2的合并成一个n的
    if(n==1) {f[0] = qpow(g[0],P-2);return;}
    get_inv(f,g,(n+1)>>1);
    int len(1);
    while(len<(n<<1)) len<<=1;
    copy(g,g+n,inv_t);
    fill(inv_t+n,inv_t+len,0);

    ntt(f,len,1);
    ntt(inv_t,len,1);
    for(int i=0;i<len;++i){
        f[i] = 1ll*f[i]*(2-1ll*f[i]*inv_t[i]%P+P)%P;
    }
    ntt(f,len,-1);
    fill(f+n,f+len,0);  // 取模
}

int f[MAXN<<2],g[MAXN<<2],dg[MAXN<<2],ans[MAXN<<2];

int main(){
    int n=read();
    for(int i=0;i<n;++i) g[i]=read();
    // 求逆
    get_inv(f,g,n);
    // 求导
    for(int i=1;i<n;++i){
        dg[i-1]=1ll*g[i]*i%P;
    }
    // 乘在一起
    int len(1);
    while(len<(n<<1)) len<<=1;
    ntt(dg,len,1),ntt(f,len,1);
    for(int i=0;i<len;++i){
        f[i] = 1ll*f[i]*dg[i]%P;
    }
    ntt(f,len,-1);
    // 再积分
    for(int i=1;i<n;++i){
        ans[i]=1ll*qpow(i,P-2)*f[i-1]%P;
    }
    for(int i=0;i<n;++i)
        printf("%d%c",ans[i],i==n-1?'\n':' ');
}

多项式指数函数

洛谷板子

法一：O(nlog^2n)分治FFT

$B(x)=e^{A(x)}$
两边求导
$B^{\prime }(x)=A^{\prime }(x)e^{A(x)}$
两边积分
$B(x)=\int A^{\prime }(x)B(x)$
发现形如分治FFT
由$A(0)=0$得$B(0)=1$（递归边界）
（1）开始的时候把A变成A’。
（2）FFT乘出来之后将a数组进行移位积分，再加到f数组中。
（3）还是积分：要在递归边界时除以i。

#include<bits/stdc++.h>

using namespace std;
#define int long long int
#define mp(a, b) make_pair(a,b)
#define vi vector<int>
#define mii map<int,int>
#define mpi map<pair<int,int>,int>
#define vp vector<pair<int,int> >
#define pb(a) push_back(a)
#define fr(i, n) for(i=0;i<n;i++)
#define rep(i, a, n) for(i0=a;i<n;i++)
#define FOR(i,a,b) for(i=(a);i<=(b);++i)
#define F first
#define S second
#define endl "\n"
#define Endl "\n"
#define fast std::ios_base::sync_with_stdio(false);cin.tie(0);cout.tie(0);
#define mod 1000000007
#define dom 998244353
#define sl(a) (int)a.length()
#define sz(a) (int)a.size()
#define all(a) a.begin(),a.end()
#define pii pair<int,int>
#define mini 2000000000000000000
#define time_taken 1.0 * clock() / CLOCKS_PER_SEC

//mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
//const long double pi = acos(-1);
//mt19937_64 mt(chrono::steady_clock::now().time_since_epoch().count());
//primes for hashing 937, 1013
template<typename T, typename U>
static inline void amin(T &x, U y) {
    if (y < x)
        x = y;
}

template<typename T, typename U>
static inline void amax(T &x, U y) {
    if (x < y)
        x = y;
}

const int MAXN = 1e5+5;
inline int read(){int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch = getchar();}while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}return x*f;}

const int P = 998244353;

inline int qpow(int x,int y){int res(1);while(y){if(y&1) res = 1ll*res*x%P;x = 1ll*x*x%P;y>>=1;}return res;}

int rev[MAXN<<2];

void change(int y[],int len){
    for(int i=0;i<len;++i){
        rev[i] = rev[i>>1]>>1;
        if(i&1) rev[i] |= len>>1;
    }
    for(int i=0;i<len;++i)
        if(i<rev[i]) swap(y[i],y[rev[i]]);
}

void ntt(int y[],int len,int on){
    change(y,len);
    for(int h=2;h<=len;h<<=1){
        int gn = qpow(3,(P-1)/h);
        for(int j=0;j<len;j+=h){
            int g = 1;
            for(int k=j;k<j+h/2;++k){
                int u = y[k];
                int t = 1ll*g*y[k+h/2]%P;
                y[k] = (u+t)%P;
                y[k+h/2] = (u-t+P)%P;
                g = 1ll*g*gn%P;
            }
        }
    }
    if(on==-1) {
        reverse(y+1,y+len);
        register int inv = qpow(len,P-2);
        for(int i=0;i<len;++i) y[i] = 1ll*y[i]*inv%P;
    }
}

int f[MAXN<<2],g[MAXN<<2];
int a[MAXN<<2],b[MAXN<<2]; // temp

void solve(int l,int r){
    if(l+1>=r) {
        // 积分：把常数除掉
        if(l)(f[l]*=qpow(l,P-2))%=P;
        else f[l]=1;
        return;
    }
    int len = (r-l)>>1; // 半长
    int mid = len+l;
    solve(l,mid); // 递归计算左区间
    memset(a+len,0,sizeof(int)*len);
    memcpy(a,f+l,sizeof(int)*len); // a中储存的是?
    memcpy(b,g,sizeof(int)*(r-l)); // b中储存的是?
    ntt(a,r-l,1),ntt(b,r-l,1);
    for(int i=0;i<r-l;++i) a[i]=a[i]*b[i]%P;
    ntt(a,r-l,-1);
    for(int i=len;i<r-l;++i) f[l+i]=(f[l+i]+a[i-1])%P;
    // 平移一项加上
    solve(mid,r); // 递归计算右区间
}


signed main() {
    fast;
    int n = read();
    for(int i=0;i<n;++i) g[i]=read();
    for(int i=0;i<n;++i) g[i-1]=g[i]*i%P;
    g[n-1]=0;
    int len(1);
    while(n>len) len<<=1;
    solve(0,len);
    for(int i=0;i<n;++i)
        printf("%d%c",f[i],i==n-1?'\n':' ');
    return 0;
}

多项式牛顿迭代

泰勒展开公式

如果$f(x)$在某个包含$x_0$的闭区间$[a,b]$中存在$n$阶导数，并在开区间$(a,b)$上存在$n+1$阶导数，则对于$[a,b]$中任意一点$x$，均存在：
$$f(x)=g(x_0)+\frac{f^{(1)}(x_0)}{1!}(x-x_0)+\frac{f^{(2)}(x_0)}{2!}(x-x_0{)}^2+\dots \dots +\frac{f^{(n)}(x_0)}{n!}(x-x_0{)}^n+R_n(x)$$
此处$R_n(x)$是一个余项。

Newton’s method

$$F(x)\equiv F_0(x)-\frac{G(F_0(x))}{G^{\prime }(F_0(x))}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$

FWT

不懂细节可以看洛谷的题解，讲的十分详细

//
// Created by artist on 2021/8/22.
//


#include<bits/stdc++.h>

using namespace std;
typedef long long ll;
#define mkp make_pair
#define fi first
#define se second
#define pb push_back
#define all(x) x.begin(),x.end()
#define DB1(args...) do { cout << #args << " : "; dbg(args); } while (0)

void dbg() { std::cout << "  #\n"; }

template<typename T, typename...Args>
void dbg(T a, Args...args) {
    std::cout << a << ' ';
    dbg(args...);
}
const int maxn = 2e5+5;
const int mod = 998244353;
ll A[maxn],B[maxn];
ll a[maxn],b[maxn];
int n;
void in(){
    for(int i=0;i<(1<<n);++i) a[i]=A[i]%mod;
    for(int i=0;i<(1<<n);++i) b[i]=B[i]%mod;
}

void get(){
    for(int i=0;i<(1<<n);++i) a[i]=a[i]*b[i]%mod;
}

void out(){
    for(int i=0;i<(1<<n);++i) printf("%lld%c",a[i],i==(1<<n)-1?'\n':' ');
}

void OR(ll f[],int x=1){
    for(int o=2,k=1;o<=(1<<n);o<<=1,k<<=1){
        for(int i=0;i<(1<<n);i+=o){
            for(int j=0;j<k;++j){
                f[i+j+k] = (f[i+j+k] + x*f[i+j])%mod;
            }
        }
    }
}

void AND(ll f[],int x=1){
    for(int o=2,k=1;o<=(1<<n);o<<=1,k<<=1){
        for(int i=0;i<(1<<n);i+=o){
            for(int j=0;j<k;++j){
                f[i+j] = (f[i+j] + x*f[i+j+k])%mod;
            }
        }
    }
}

void XOR(ll f[],int x=1){
    for(int o=2,k=1;o<=(1<<n);o<<=1,k<<=1){
        for(int i=0;i<(1<<n);i+=o){
            for(int j=0;j<k;++j){
                f[i+j] = (f[i+j] + f[i+j+k])%mod;
                f[i+j+k] = (f[i+j] - f[i+j+k]+mod - f[i+j+k]+mod)%mod;
                f[i+j]=f[i+j]*x%mod,f[i+j+k]=f[i+j+k]*x%mod;
            }
        }
    }
}

signed main() {
    scanf("%d",&n);
    for(int i=0;i<(1<<n);++i) scanf("%lld",&A[i]);
    for(int i=0;i<(1<<n);++i) scanf("%lld",&B[i]);
    in();OR(a,1);OR(b,1);get();OR(a,mod-1);out();
    in();AND(a,1);AND(b,1);get();AND(a,mod-1);out();
    in();XOR(a,1);XOR(b,1);get();XOR(a,(mod+1)>>1);out();
}

生成函数

封闭形式和形式幂级数形式之间的转化。

普通生成函数

$F(x)={\sum }_{n\ge 0}\left(\genfrac{}{}{0ex}{}{m+n}{n}\right)x^n=\frac{1}{{(1-x)}^{m+1}}$。
证明：归纳法（见OI wiki）
有一种常见的转化操作：（变换贡献指向）
例子1：
${\sum }_{n\ge 0}{\sum }_{i=0}^n\left(\genfrac{}{}{0ex}{}{n}{i}\right)x^{n+i}$
$={\sum }_{n\ge 0}{x}^n{\sum }_{i=0}^n\left(\genfrac{}{}{0ex}{}{n-i}{i}\right)$
例子2：
${\sum }_{n\ge 0}\left(\genfrac{}{}{0ex}{}{m+n-1}{n}\right)x^n{\sum }_{n\ge 0}{x}^n$
$={\sum }_{n\ge 0}{x}^n{\sum }_{i=0}^n\left(\genfrac{}{}{0ex}{}{m+i-1}{i}\right)$
例子3：（卡特兰数）

中的第三到第四步：
$=1+{\sum }_{n\ge 1}{\sum }_{i=0}^{n-1}{H}_i{H}_{n-i-1}{x}^n$
$=1+{\sum }_{n\ge 1}{x}^n{\sum }_{i=0}^{n-1}{H}_i{H}_{n-i-1}$
从 贡献给x 变成 从x贡献给x+i+1，同时转化为卡特兰数的生成函数的形式。
$=1+x{\sum }_{i\ge 0}{H}_i{x}^i{\sum }_{n\ge 0}{H}_n{x}_n$

分数多项式封闭形式通法：待定系数法

多项式快速幂板子

//
// Created by Artist on 2021/10/23.
//

#include<bits/stdc++.h>

using namespace std;
typedef long long ll;
#define mkp make_pair
#define fi first
#define se second
#define pb push_back
#define all(x) x.begin(),x.end()
#define DB1(args...) do { cout << #args << " : "; dbg(args); } while (0)
#define endl '\n'
#define reg register

void dbg() { std::cout << "  #\n"; }

template<typename T, typename...Args>
void dbg(T a, Args...args) {
    std::cout << a << ' ';
    dbg(args...);
}

void io() { ios::sync_with_stdio(0), cin.tie(0), cout.tie(0); }
const int maxn = 1e5+5;
const int P = 998244353;
const int g = 3;
string s;

inline int read() {
    reg char ch=getchar();
    reg int x=0,f=1;
    while(ch<'0'||ch>'9') {
        if(ch=='-') f=-1;
        ch=getchar();
    }
    while(isdigit(ch)) {
        x=(x<<3)+(x<<1)+ch-'0';
        ch=getchar();
    }
    return x*f;
}

inline int readm() {
    reg char ch=getchar();
    reg ll x=0,f=1;
    while(ch<'0'||ch>'9') {
        if(ch=='-') f=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9') {
        x=((x<<3)+(x<<1)+ch-'0')%P;
        ch=getchar();
    }
    return x*f;
}


inline int qpow(int a,int n) {
    int ans=1;
    while(n) {
        if(n&1) ans=1ll*ans*a%P;
        a=1ll*a*a%P;
        n>>=1;
    }
    return ans;
}

void ntt(int a[],int n,int t) {
    static int rev[maxn<<2];
    for(int i=1;i<n;++i) {
        rev[i] = (rev[i>>1]>>1)|(i&1?n>>1:0);
        if(i>rev[i]) swap(a[i],a[rev[i]]);
    }
    for(int i=2;i<=n;i<<=1) {
        int wn=qpow(g,(P-1)/i);
        if(t==-1) wn=qpow(wn,P-2);
        for(int j=0;j<n;j+=i) {
            int w=1;
            for(int k=j;k<j+i/2;++k) {
                int u=a[k];
                int v=1ll*a[k+i/2]*w%P;
                a[k]=(1ll*u+v)%P;
                a[k+i/2]=(1ll*u-v+P)%P;
                w=1ll*w*wn%P;
            }
        }
    }
    if(t==-1) {
        ll inv=qpow(n,P-2);
        for(int i=0;i<n;++i) a[i]=1ll*a[i]*inv%P;
    }
}

void polyinv(int f[],int h[],int n) {
    static int tmp[maxn<<2];
    int len;
    for(len=1;len<n;len<<=1);
    fill(f,f+len+len+2,0);
    fill(tmp,tmp+len+1,0);
    f[0] = qpow(h[0],P-2);
    for(int t=2;t<=len;t<<=1) {
        const int t2=t<<1;
        copy(h,h+t,tmp);
        fill(tmp+t,tmp+t2,0);
        ntt(f,t2,1);
        ntt(tmp,t2,1);
        for(int i=0;i!=t2;++i)
            f[i]=1ll*f[i]*(2-1ll*f[i]*tmp[i]%P+P)%P;
        ntt(f,t2,-1);
        fill(f+t,f+t2,0);
    }
    fill(f+n,f+len+1,0);
}

void polyln(int f[],int h[],int n) {
    static int tmp2[maxn<<2],tmp3[maxn<<2];
    int lim;
    for(lim=1;lim<(n<<1);lim<<=1);
    fill(tmp2,tmp2+lim,0);
    fill(tmp3,tmp3+lim,0);
    fill(f,f+lim,0);
    polyinv(tmp3,h,n);
    for(int i=1;i<n;++i) tmp2[i-1]=1ll*h[i]*i%P;
    ntt(tmp2,lim,1);
    ntt(tmp3,lim,1);
    for(int i=0;i<lim;++i) tmp3[i]=1ll*tmp3[i]*tmp2[i]%P;
    ntt(tmp3,lim,-1);
    for(int i=1;i<n;++i) f[i]=1ll*qpow(i,P-2)*tmp3[i-1]%P;
    fill(f+n,f+lim,0);
}

void polyexp(int f[],int h[],int n) {
    static int tmp4[maxn<<2];
    if(n==1) {f[0]=1;return;}
    int lim;
    for(lim=1;lim<(n<<1);lim<<=1);
    fill(tmp4,tmp4+lim,0);
    polyexp(f,h,n+1>>1),polyln(tmp4,f,n);
    tmp4[0]=(1ll*h[0]+1-tmp4[0]+P)%P;
    for(int i=1;i<n;++i) tmp4[i]=(1ll*h[i]-tmp4[i]+P)%P;
    ntt(tmp4,lim,1),ntt(f,lim,1);
    for(int i=0;i<lim;++i) f[i]=1ll*f[i]*tmp4[i]%P;
    ntt(f,lim,-1);
    fill(f+n,f+lim,0);
}

int A[maxn<<2],B[maxn<<2],C[maxn<<2];
signed main() {
    int n=read(),k=readm();
    for(int i=0;i<n;++i) A[i]=read();
    polyln(B,A,n);
    for(int i=0;i<n;++i) B[i]=1ll*B[i]*k%P;
    polyexp(C,B,n);
    for(int i=0;i<n;++i) printf("%d ",C[i]);
}