模意义下的FFT算法
//写在前面 单就FFT算法来说的话,下面只给出个人认为比较重要的推导,详细的介绍可参考 FFT算法学习笔记
令v[n]是长度为2N的实序列,V[k]表示该实序列的2N点DFT。定义两个长度为N的实序列g[n]和h[n]为
g[n]=v[2n], h[n]=v[2n+1], 0<=n<N
则可进行如下推导
这里所用的FFT算法能够实现O(nlogn)复杂度的离散傅里叶变换和上面最后所得的关系密切相关。
下面进入正题——模意义下的FFT
还是需要先了解一下关于 复序列的DFT的对称性质及一些补充定义
由此,可以试想,假设说要模的素数p为1e8级别大小,那么我们可以把原始的实序列x[n]“拆”一下。
下面假设我们要求的是x[n]卷积y[n]的结果t[n]。
假设q是sqrt(p)级别的一个数,我们可以把x[n]/q存到复序列x1[n]的实部,x[n]%q存到复序列x1[n]的虚部。这时,对x1[n]、y1[n]求DFT,再由X1[k]*Y1[k]得到T1[k],整个运算过程中能够产生的最大浮点数为N*q^2级别,一般来说还是在可以接受的范围内的。
接下来考虑从卷积结果{T1[k]}中恢复出原始的t[n]的过程。
看一下T1[k]的组成
到这里差不多就可以结束了。发现上面最后一行等号右边有四个“乘积”,我们可以把上面四个乘积分别单独拿出来,求IDFT就可以恢复出x/y_re/im卷积的结果,之后针对不同的乘积,考虑需要在模p意义下乘上q^2、q^1或q^0,来进行恢复就可以了。
奉上模板
namespace FFT_MO //前面需要有 mod(1e8~1e9级别),upmo(a,b) 的定义 {const int FFT_MAXN=1<<18;const db pi=3.14159265358979323846264338327950288L;struct cp{db a,b;cp(double a_=0,double b_=0){a=a_,b=b_;}cp operator +(const cp&rhs)const{return cp(a+rhs.a,b+rhs.b);}cp operator -(const cp&rhs)const{return cp(a-rhs.a,b-rhs.b);}cp operator *(const cp&rhs)const{return cp(a*rhs.a-b*rhs.b,a*rhs.b+b*rhs.a);}cp operator !()const{return cp(a,-b);}}nw[FFT_MAXN+1],f[FFT_MAXN],g[FFT_MAXN],t[FFT_MAXN]; //a<->f,b<->g,t<~>c int bitrev[FFT_MAXN]; void fft_init() //初始化 nw[],bitrev[] {int L=0;while((1<<L)!=FFT_MAXN) L++;for(int i=1;i<FFT_MAXN;i++) bitrev[i]=bitrev[i>>1]>>1|((i&1)<<(L-1));for(int i=0;i<=FFT_MAXN;i++) nw[i]=cp((db)cosl(2*pi/FFT_MAXN*i),(db)sinl(2*pi/FFT_MAXN*i));}// n已保证是2的整数次幂 // flag=1:DFT | flag=-1: IDFTvoid dft(cp *a,int n,int flag=1) {int d=0;while((1<<d)*n!=FFT_MAXN) d++;for(int i=0;i<n;i++) if(i<(bitrev[i]>>d))swap(a[i],a[bitrev[i]>>d]);for(int l=2;l<=n;l<<=1){int del=FFT_MAXN/l*flag; // 决定 wn是在复平面是顺时针还是逆时针变化,以及变化间距 for(int i=0;i<n;i+=l){cp *le=a+i,*ri=a+i+(l>>1);cp *w=flag==1? nw:nw+FFT_MAXN; // 确定wn的起点 for(int k=0;k<(l>>1);k++){cp ne=*ri * *w;*ri=*le-ne,*le=*le+ne;le++,ri++,w+=del;}}}if(flag!=1) for(int i=0;i<n;i++) a[i].a/=n,a[i].b/=n;}// convo(a,n,b,m,c) a[0..n]*b[0..m] -> c[0..n+m]void convo(LL *a,int n,LL *b,int m,LL *c){for(int i=0;i<=n+m;i++) c[i]=0;int N=2;while(N<=n+m) N<<=1; // N是c扩展后的长度 for(int i=0;i<N;i++) //扩展 a[],b[] ,存入f[],g[],注意取模 {LL aa=i<=n?a[i]:0,bb=i<=m? b[i]:0; aa%=mod,bb%=mod;f[i]=cp(db(aa>>15),db(aa&32767));g[i]=cp(db(bb>>15),db(bb&32767));}dft(f,N),dft(g,N);for(int i=0;i<N;i++) // 恢复虚部两个“乘积”(乘积具体意义见上文) {int j=i? N-i:0;t[i]=((f[i]+!f[j])*(!g[j]-g[i])+(!f[j]-f[i])*(g[i]+!g[j]))*cp(0,0.25);}dft(t,N,-1);for(int i=0;i<=n+m;i++) upmo(c[i],(LL(t[i].a+0.5))%mod<<15);for(int i=0;i<N;i++) // 恢复实部两个“乘积” {int j=i? N-i:0;t[i]=(!f[j]-f[i])*(!g[j]-g[i])*cp(-0.25,0)+cp(0,0.25)*(f[i]+!f[j])*(g[i]+!g[j]);}dft(t,N,-1);for(int i=0;i<=n+m;i++) upmo(c[i],LL(t[i].a+0.5)+(LL(t[i].b+0.5)%mod<<30));} }
模板
举个栗子~ hdu 6088 Rikka with Rock-paper-scissors (2017 多校第五场 1004) 【组合数学 + 数论 + 模意义下的FFT】
//本博客主要参考资料:数字信号处理——基于计算机的方法(第四版) [美] Sanjit K. Mitra 著 余翔宇 译
转载请注明出处 http://www.cnblogs.com/Just--Do--It/p/7892254.html
谢谢阅读
转载于:https://www.cnblogs.com/Just--Do--It/p/7892254.html
模意义下的FFT算法相关推荐
- 10.24T3 解方程 取模意义下运算+秦九韶算法
#1228 解方程 描述 已知多项式方程: a0+a1x+a2x^2+..+anx^n=0 求这个方程在[1, m ] 内的整数解(n 和m 均为正整数) 输入 输入共n + 2 行. 第一行包含2 ...
- CERC2015 Frightful Formula 神奇的模意义下分数
上网查了下这道题的正解是FFT......然而机智的xushu用待定系数法A了.... f[i,j]+k=a(f[i][j-1]+k)+b(f[i-1,j]+k) 解得 k=c/(a+b-1) 于是令 ...
- Codeforces 1106F Lunar New Year and a Recursive Sequence 矩阵快速幂,原根转化模意义下对数,BSGS
文章目录 题意 题解 对数法转指数线性递推 原根与模意义下求对数 拔山盖世! 最终步骤 Problem Origin 狠搞了一个多星期,做出来之后仍然一知半解,写个博客重理思路. 题意 定义序列fff ...
- Newcoder Wannafly13 B Jxy军训(费马小定理、分数在模意义下的值)
链接:https://www.nowcoder.com/acm/contest/80/B 题目描述 在文某路学车中学高一新生军训中,Jxc正站在太阳下站着军姿,对于这样的酷热的阳光,Jxc 表示非常不 ...
- 求解斐波那契数列模$p$意义下最短循环节
如题,毕克老师给我们出的noip(NOIplus)模拟赛的\(Day1T1\) 首先我们知道斐波那契数列的特征根 \[\phi_1=\frac{1+\sqrt{5}}{2}\] \[\phi_2=\f ...
- 算法复杂性和渐进复杂意义下的记号 O、Ω、θ
算法复杂性 算法复杂性 === 算法运行时所需要的计算机资源的量 通常指的是空间.时间资源 影响时间复杂性的因素 影响时间复杂性的因素有: 问题规模nnn.输入序列III.算法本身AAA 问题规模:如 ...
- 深入浅出理解FFT算法。通俗易懂,xilinxIP核仿真
深入浅出理解FFT算法,通俗易懂,用xilinxIP核心仿真 1.前言:傅里叶变换:时域到频域的转换 FS连续时间周期傅里叶级数->DFS离散傅里叶级数->FT连续时间非周期信号的傅里叶变 ...
- FFT算法8点12位硬件实现 (verilog)
FFT算法8点12位硬件实现 (verilog) 1 一.功能描述: 1 二.设计结构: 2 三.设计模块介绍 3 1.蝶形运算(第一级) 3 2.矢量角度旋转(W) 4 3.CORDIC 结果处理 ...
- FFT算法的DSP实现
FFT算法移植到DSP的过程 一, 学习导向 为什么要学习5000系列DSP?我想这就像为什么要学习<信号与系统>,<模拟电路>一样,他不仅是理论性强的东西,更是直接具 ...
最新文章
- Android Build.VERSION.SDK_INT
- jquery radio 取值
- R语言message函数、warning()函数和stop()函数输出程序运行健康状态信息实战
- IP BASE对OSPF的支持版本
- MOSS007 服务器的配置
- 上一局APP玩边画边猜,第1次见人使用道具,我的游戏体验上升了
- SAP UI5 Label related stuff and accessibility研究
- 使用org.apache.commons.io.FileUtils,IOUtils工具类操作文件
- Linux命令之 chsh -- 用来更换登录系统时使用的shell
- CriminalIntent项目开发
- 单例模式中的线程安全问题
- table 样式美化
- Java使用Spire.pdf提取PDF中想要的图片
- Hdu4939 Stupid Tower Defense
- 猴子意念打字,有可能敲出莎士比亚全集
- 文笔很差系列1 - 也谈谈AlphaGo
- 27年,微软IE结束了!
- Alibaba与gofair的对比
- Qt入门开发__贪吃蛇小游戏
- 流畅的python学习笔记
热门文章
- c语言变量ppt,C语言程序设计-变量.ppt
- php冒泡和选择排序,选择排序vs冒泡排序
- dingo php,详细介绍Laravel+Dingo/Api 自定义响应
- mysql数据库导入导出
- php dirtoarray,PHP Ds\Stack toArray()用法及代码示例
- 研发项目管理中需注意的人性弱点(Z)
- vasp 5.2编译方法
- 随机森林分类器_建立您的第一个随机森林分类器
- 在什么情况下,刘强东会丧失京东的控制权?
- 单片机小白学步系列(十七) 单片机/计算机系统概述:核心模块