分治FFT

最新推荐文章于 2021-08-25 18:11:14 发布

weixin_30955341

最新推荐文章于 2021-08-25 18:11:14 发布

阅读量127

点赞数

原文链接：http://www.cnblogs.com/guoshaoyang/p/11282228.html

版权

分治FFT

概述

严格的分治FFT为“一般形式”，而我们把所有带有FFT的分治都称为分治FFT，所以这个名称并没有什么意义。
分治FFT代码比较复杂，低阶的比赛里应该不会出现。

一般形式

严格的分治FFT负责解决类似\(\displaystyle b_k=\sum^{k-1}_{i=1}{b_ia_{k-i}}\)的数列，其中整个式子的形状符合卷积，且所求数列的每一项都依赖于前面的项。

已知数组\(a\)，求数组\(b\)，其中\(\displaystyle b_k=\sum^{k-1}_{i=1}{b_ia_{k-i}}\)

考虑分治，对于区间\([l,r)\)，先计算\([l,m)\)的答案。对于\(k\geq m\)，将\(b_k\)拆成\(\displaystyle \sum^{m-1}_{i=0}b_ia_{k-i}\)和\(\displaystyle \sum^{k-1}_{i=m}b_ia_{k-i}\)两部分，后者可以给数组下标添加偏移量之后递归计算，而前者可以通过普通的卷积计算。

总时间复杂度为\(\Theta(n \log^2 n)\)。

实现

假设我们求出了\(l\to mid\)的答案，要求这些点对\(mid+1\to r\)的影响，那么对右半边点\(x\)的贡献为:\(\displaystyle w_x=\sum_{i=l}^{mid} f[i] * g[x-i]\)这部分可以利用卷积来快速计算。计算完以后，答案直接加到答案数组就可以了。需要注意的是，如果要求左边点对右边点的影响，首先整个区间以左对该区间的贡献应该先求出。所以分治过程为先分治左边，在求出中间，然后在递归右边。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int MAXN=1e6+10;
const LL MOD=998244353,G=3,iG=332748118;
inline LL fpm(LL base,LL p){
    LL ret=1;
    while(p){
        if(p&1)
            ret=ret*base%MOD;
        base=base*base%MOD;
        p>>=1;
    }
    return ret;
}
void exgcd(LL a,LL b,LL &x,LL &y){
    if(!b){
        x=1;
        y=0;
        return;
    }
    exgcd(b,a%b,y,x);
    y-=x*(a/b);
}
inline LL inv(LL x){
    LL ii,jj;
    exgcd(x,MOD,ii,jj);
    return (ii%MOD+MOD)%MOD;
}
int lim,rev[MAXN],N,L;
inline void prec(int l,int r){
    lim=1;
    L=0;
    while(lim<r-l){
        L++;
        lim<<=1;
    }
    for(int i=0;i<lim;i++)
        rev[i]=rev[i>>1]>>1|(i&1)<<(L-1);
}
inline void NTT(LL *a,int type){
    for(int i=0;i<lim;i++)
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    for(int hf=1;hf<lim;hf<<=1){
        int len=hf<<1;
        LL Wn=fpm(type==1?G:iG,(MOD-1)/len);
        for(int j=0;j<lim;j+=len){
            LL w=1;
            for(int k=0;k<hf;k++){
                LL t1=a[j+k],t2=a[j+k+hf]*w%MOD;
                a[j+k]=(t1+t2)%MOD;
                a[j+k+hf]=(t1-t2)%MOD;
                w=w*Wn%MOD;
            }
        }
    }
}
LL f[MAXN],g[MAXN],A[MAXN],B[MAXN];
inline void cdq(int l,int r){
    if(l+1>=r||l>=N)
        return;
    int mid=(l+r)>>1;
    cdq(l,mid);
    prec(l,r);
    memcpy(A,f+l,((r-l)/2)<<3);
    memcpy(B,g,(r-l)<<3);
    memset(A+(r-l)/2,0,((r-l)/2)<<3);
    NTT(A,1);
    NTT(B,1);
    for(int i=0;i<r-l;i++)
        A[i]=A[i]*B[i]%MOD;
    NTT(A,-1);
    int ii=inv(r-l);
    for(int i=0;i<r-l;i++)
        A[i]=A[i]*ii%MOD;
    for(int i=(r-l)>>1;i<r-l;i++)
        f[i+l]=(f[i+l]+A[i])%MOD;
    cdq(mid,r);
}
int main(){
    scanf("%d",&N);
    for(int i=1;i<N;i++)
        scanf("%lld",g+i);
    prec(0,N);
    f[0]=1;
    cdq(0,lim);
    for(int i=0;i<N;i++)
        printf("%lld ",(f[i]+MOD)%MOD);
    return 0;
}

例题

P4721 【模板】分治 FFT

见“实现”

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int MAXN=1e6+10;
const LL MOD=998244353,G=3,iG=332748118;
inline LL fpm(LL base,LL p){
    LL ret=1;
    while(p){
        if(p&1)
            ret=ret*base%MOD;
        base=base*base%MOD;
        p>>=1;
    }
    return ret;
}
void exgcd(LL a,LL b,LL &x,LL &y){
    if(!b){
        x=1;
        y=0;
        return;
    }
    exgcd(b,a%b,y,x);
    y-=x*(a/b);
}
inline LL inv(LL x){
    LL ii,jj;
    exgcd(x,MOD,ii,jj);
    return (ii%MOD+MOD)%MOD;
}
int lim,rev[MAXN],N,L;
inline void prec(int l,int r){
    lim=1;
    L=0;
    while(lim<r-l){
        L++;
        lim<<=1;
    }
    for(int i=0;i<lim;i++)
        rev[i]=rev[i>>1]>>1|(i&1)<<(L-1);
}
inline void NTT(LL *a,int type){
    for(int i=0;i<lim;i++)
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    for(int hf=1;hf<lim;hf<<=1){
        int len=hf<<1;
        LL Wn=fpm(type==1?G:iG,(MOD-1)/len);
        for(int j=0;j<lim;j+=len){
            LL w=1;
            for(int k=0;k<hf;k++){
                LL t1=a[j+k],t2=a[j+k+hf]*w%MOD;
                a[j+k]=(t1+t2)%MOD;
                a[j+k+hf]=(t1-t2)%MOD;
                w=w*Wn%MOD;
            }
        }
    }
}
LL f[MAXN],g[MAXN],A[MAXN],B[MAXN];
inline void cdq(int l,int r){
    if(l+1>=r||l>=N)
        return;
    int mid=(l+r)>>1;
    cdq(l,mid);
    prec(l,r);
    memcpy(A,f+l,((r-l)/2)<<3);
    memcpy(B,g,(r-l)<<3);
    memset(A+(r-l)/2,0,((r-l)/2)<<3);
    NTT(A,1);
    NTT(B,1);
    for(int i=0;i<r-l;i++)
        A[i]=A[i]*B[i]%MOD;
    NTT(A,-1);
    int ii=inv(r-l);
    for(int i=0;i<r-l;i++)
        A[i]=A[i]*ii%MOD;
    for(int i=(r-l)>>1;i<r-l;i++)
        f[i+l]=(f[i+l]+A[i])%MOD;
    cdq(mid,r);
}
int main(){
    scanf("%d",&N);
    for(int i=1;i<N;i++)
        scanf("%lld",g+i);
    prec(0,N);
    f[0]=1;
    cdq(0,lim);
    for(int i=0;i<N;i++)
        printf("%lld ",(f[i]+MOD)%MOD);
    return 0;
}

转载于:https://www.cnblogs.com/guoshaoyang/p/11282228.html

weixin_30955341

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
分治FFT

分治FFT概述严格的分治FFT为“一般形式”，而我们把所有带有FFT的分治都称为分治FFT，所以这个名称并没有什么意义。分治FFT代码比较复杂，低阶的比赛里应该不会出现。一般形式严格的分治FFT负责解决类似\(\displaystyle b_k=\sum^{k-1}_{i=1}{b_ia_{k-i}}\)的数列，其中整个式子的形状符合卷积，且所求数列的每一项都依赖于前面的项。...
复制链接

扫一扫