首页 > 其他分享 >分治FFT

分治FFT

时间:2023-07-13 21:14:41浏览次数:30  
标签:FFT int res rep 分治 mid inline mod

考虑一下递推式:

\[ f_{i}=\sum\limits_{j=1}^i g_jf_{i-j} \]

如果要暴力计算的话复杂度是 \(n^2\),考虑到后面的是卷积的形式,但是唯一的问题是 \(f\) 自己得出自己。考虑分治。设当前分治区间为 \(l,r\),首先分治计算 \(l,mid\) 的所有 \(f\) 值,然后考虑 \(l,mid\) 给 \(mid+1,r\) 的 \(f\) 值的贡献,具体来说,我们把 \(l,mid\) 的 \(f\) 和 \(0,r-l\) 的 \(g\) 卷一下得到给后面的贡献。

考虑这个题也可以求逆来做,设 \(F(x)=\sum f_ix^i\),设 \(G(x)=\sum g_ix^i\)。由于 \(f_0=1\),\(g_0=0\),考虑有 \(F(x)=F(x)G(x)+f_0\Rightarrow F(x)\equiv \frac{f_0}{1-G(x)} \bmod x^n\)。求逆即可。

分治 NTT:

int n,g[N],tr[N],F[N],G[N],H[N],f[N];

inline int ksm(int a,int b,int mod){int res=1;while(b){if(b&1)res=1ll*a*res%mod;a=1ll*a*a%mod;b>>=1;}return res;}
inline int inv(int a){return ksm(a,mod-2,mod);}
inline void Gettr(int n){
    for(int i=0;i<n;i++) tr[i]=(tr[i>>1]>>1)|((i&1)?(n>>1):0);
}
inline void NTT(int *f,int n,int op){
    rep(i,0,n-1) if(i<tr[i]) swap(f[i],f[tr[i]]);
    for(int i=2;i<=n;i<<=1){
        int x=ksm(gg,(mod-1)/i,mod),l=i>>1;if(op==0) x=inv(x);
        for(int j=0;j<n;j+=i){
            int now=1;
            for(int k=j;k<j+l;k++){
                int tt=1ll*now*f[k+l]%mod;
                f[k+l]=(f[k]-tt)%mod;f[k]=(f[k]+tt)%mod;
                now=1ll*now*x%mod;
            }
        }
    }
    if(op==0){
        int invn=inv(n);rep(i,0,n-1) f[i]=1ll*f[i]*invn%mod;
    }
}
inline void Solve(int l,int r){
    if(l==r){if(l==0)f[0]=1;return;}
    int mid=(l+r)>>1,len=r-l+1;
    Solve(l,mid);rep(i,0,len-1) F[i]=G[i]=0;
    rep(i,0,mid-l) F[i]=f[i+l];
    rep(i,0,len-1) G[i]=g[i];
    // printf("l=%d r=%d\n",l,r);
    // rep(i,0,len-1) printf("F[%d]=%d\n",i,F[i]);
    // rep(i,0,len-1) printf("G[%d]=%d\n",i,G[i]);
    int nl=1;while(nl<len*2){nl<<=1;}Gettr(nl);
    rep(i,mid-l+1,nl-1) F[i]=0;rep(i,len,nl-1) G[i]=0;
    NTT(F,nl,1);NTT(G,nl,1);
    // rep(i,0,len-1) printf("F[%d]=%d\n",i,F[i]);
    // rep(i,0,len-1) printf("G[%d]=%d\n",i,G[i]);
    rep(i,0,nl-1) F[i]=1ll*F[i]*G[i]%mod;
    NTT(F,nl,0);rep(i,mid+1-l,r-l) (f[i+l]+=F[i])%=mod;
    // rep(i,mid+1,r-l) printf("F[%d]=%d\n",i,F[i]);
    Solve(mid+1,r);
}

int main(){
    // freopen("my.in","r",stdin);
    // freopen("my.out","w",stdout);
    read(n);rep(i,1,n-1) read(g[i]);
    Solve(0,n-1);
    rep(i,0,n-1) printf("%d ",(f[i]+mod)%mod);
    return 0;
}

求逆:

int tr[N];
inline int ksm(int a,int b,int mod){int res=1;while(b){if(b&1) res=1ll*a*res%mod;a=1ll*a*a%mod;b>>=1;}return res;}

inline void Gettr(int n){
    for(int i=0;i<n;i++) tr[i]=(tr[i>>1]>>1)|((i&1)?(n>>1):0);
}
inline void NTT(int *f,int len,int flag){
    for(int i=0;i<len;i++) if(i<tr[i]) swap(f[i],f[tr[i]]);
    for(int p=2;p<=len;p<<=1){
        int md=ksm(g,(mod-1)/p,mod),l=p>>1;
        if(flag==-1) md=ksm(md,mod-2,mod);
        for(int i=0;i<len;i+=p){
            int buf=1;
            for(int j=i;j<i+l;j++){
                int tt=1ll*f[j+l]*buf%mod;
                f[j+l]=((f[j]-tt)%mod+mod)%mod;
                f[j]=(f[j]+tt)%mod;buf=1ll*buf*md%mod;
            }
        }
    }
}

int n,m,c[N],a[N],b[N];

inline void GetInv(int len,int *a,int *b){
    if(len==1){b[0]=ksm(a[0],mod-2,mod);return;}
    GetInv((len+1)>>1,a,b);m=1;while(m<(len<<1)) m<<=1;
    Gettr(m);for(int i=0;i<len;i++) c[i]=a[i];
    for(int i=len;i<m;i++) c[i]=0;
    // printf("m=%d\n",m);
    // for(int i=0;i<m;i++) printf("%d ",c[i]);
    NTT(c,m,1);NTT(b,m,1);
    for(int i=0;i<m;i++) b[i]=1ll*(2-1ll*b[i]*c[i]%mod+mod)%mod*b[i]%mod;
    NTT(b,m,-1);int inv=ksm(m,mod-2,mod);for(int i=0;i<m;i++) b[i]=1ll*b[i]*inv%mod;
    for(int i=len;i<m;i++) b[i]=0;
    // printf("len=%d\n",len);
    // for(int i=0;i<n;i++) printf("%d ",b[i]);puts("");
}

int main(){
    // freopen("my.in","r",stdin);
    // freopen("my.out","w",stdout);
    read(n);for(int i=1;i<n;i++) read(a[i]);
    rep(i,0,n-1) a[i]=-a[i];a[0]++;
    GetInv(n,a,b);
    rep(i,0,n-1) printf("%d ",b[i]);
    return 0;
}

标签:FFT,int,res,rep,分治,mid,inline,mod
From: https://www.cnblogs.com/HeNuclearReactor/p/17552189.html

相关文章

  • CDQ 分治
    CDQ分治的思想最早由IOI2008金牌得主陈丹琦在高中时整理并总结,因此得名。适用范围多用于统计有特征的点对\((i,j)\)的数量或找出有特征的点对。流程对于一个区间\([L,R]\)中的点对\((i,j)\)(\(L\lei<j\leR\)),考虑分为三部分(\(mid=\frac{L+R}2\)):\(L\lei<j\lemid\)......
  • cdq分治学习笔记
    做着做着cdq分治的题发现自己太菜了写不出来题XD,所以来写写学习笔记。CDQ分治CDQ分治一般有两个用途:解决偏序问题、将动态问题转化为静态问题。偏序问题我们先来介绍第一个问题:偏序问题。普通的二维偏序就类似于逆序对,每个元素有两个属性\(a_i,b_i\),我们需要统计\(a_......
  • matlab傅里叶变换FFT,自编的fft对不足位进行补0, 频谱图和相位图去下,已对幅值进行修正
    matlab傅里叶变换FFT,自编的fft对不足位进行补0,频谱图和相位图去下,已对幅值进行修正。ID:6925626214503643......
  • 6307: 网线主管 二分/分治
    描述 仙境的居民们决定举办一场程序设计区域赛。裁判委员会完全由自愿组成,他们承诺要组织一次史上最公正的比赛。他们决定将选手的电脑用星形拓扑结构连接在一起,即将它们全部连到一个单一的中心服务器。为了组织这个完全公正的比赛,裁判委员会主席提出要将所有选手的电脑等距离......
  • 【树分治】HDOJ 5016
    先预处理出每个点的最近点。。。。然后考虑固定根,对于每个根,从根到子树的路径中任选两条,如果dis[i]+dis[j]<w[j]。那么i就是j的一个合法点。。。。然后递归处理子树即可。。。。扩栈才过的。。。#include<iostream>#include<queue>#include<stack>#include<map>#includ......
  • 线段树分治 学习笔记
    离线算法。在时间轴上建线段树(可能要事先离散化),要维护的东西用vector什么的挂在线段树的节点上,DFS一遍线段树,每次进入一个节点就加入要维护的东西,离开时撤销即可。由于DFS的特性,只需支持最近的undo,用stack可维护。......
  • 快速傅里叶变换(FFT)学习笔记
    有关多项式的一个基础算法,学起来比较困难。快速傅里叶变换和傅里叶变换没什么关系,也不是傅里叶发明的。这种算法用于在\(O(n\logn)\)时间复杂度内求出两个多项式的卷积(相当于多项式相乘)。前置知识多项式的表示\(n\)项式等价于\(n-1\)次项式。(每个次项的系数都不为零)系......
  • CDQ分治 学习笔记
    按@ouuan大佬所说,CDQ分治可以当作ex归并看待。它的思想和归并排序十分相似:假设要对区间\([l,r)\)处理先不管\([\text{mid},r)\),计算\([l,mid)\)同理计算\([mid,r)\)补回之前忽略的部分,即“归并”例:三维偏序给定\(n\)个点\((a,b,c)\),求\(a_1\lea_2\we......
  • Matlab对wav文件做fft分析
    1.代码%指定要读取的.wav文件路径filename='jay.wav';%使用audioread函数读取.wav文件[sound_data,sample_rate]=audioread(filename);sound_data=sound_data(:,1);%计算音频数据的长度sound_length=length(sound_data);%计算FFT的点数%fft_points=......
  • 分治专题
    在牛子老师的博客下边看到yspm给了CF1019E。看了一眼,不会。看了题解,我超边分治+闵可夫斯基和,一个都不会。乐。还有20天,还能补多少坑呢,不好说。仍然是每天高压作业。但是出乎意料的晚上不是很失眠,虽然说醒了以后还是很困。现象:让大象出现的事物或者方法。大象是一种体量......