FFT入门

news/2025/2/4 7:41:34/文章来源:https://www.cnblogs.com/hzoi-Cu/p/18697705

前置知识

复数 也可以参考高中数学课本,这里只会介绍 fft 需要的(默认已经入门复数)。

多项式的相关概念。

点值表示法:假设 \(f(x)\) 是一个 \(n-1\) 次多项式,那么将 \(n\)不同的 \(x\) 代入,可以得到 \(n\)\(y\)。这 \(n\) 个点对 \((x,y)\) 唯一确定了该多项式。那么就可以通过多项式求出其点值表示,也可以反过来。

注意:以下 \(n\) 均为 \(2\) 的整数次幂,若不足则补零(显然不会影响结果)

FFT 简介

用来加速多项式乘法的一个东西。而普通多项式乘法是 \(O(n^2)\) 的。但是如果是点值表示法,则是 \(O(n)\) 的(比如 \(c(x) = a(x)\times b(x)\),那么只需要枚举 \(2\times n\) 个不同的 \(x\) 即可求出 \(c\) 的点值表示)。

考虑如何快速将两个多项式 \(a(x),b(x)\) 转成点值表示,再将一个多项式 \(c(x)\) 从点值表示转化回来即可。

这里就用到 FFT 了。

离散傅里叶变换

其实就是朴素版 FFT。

就是上面代入的 \(n\)\(x\)\(n\) 个复数。但这 \(n\) 个复数不是随便找的,而是 \(n\) 次单位根。

简单介绍一下,就是 \(x^n=1\) 在复数意义下所有的根。显然这玩意有 \(n\) 个。将这几个根按幅角从小到大排序,\(0\) 开始编号,第 \(i\)\(n\) 次单位根记为 \(w_n^i\)。没有特殊声明时,一般特指第一个单位根,简记为 \(w_n\)

可以算出 \(w_n^k=\cos(\frac{2k\pi}{n})+i\sin(\frac{2k\pi}{n})\)

有两个性质:

  1. \(w_{2n}^{2k}=w_n^k\)。画个图理解一下。
  2. \(w_{n}^{k+\frac{n}{2}}-w_{n}^{k}\),画图发现它们关于原点对称。

画个图出来,大概这样子。

image

好看是好看, 但是为什么非要选择这 \(n\) 个点呢?

有个结论,将 \(a(x)\) 的离散傅里叶变换的结果作为 \(b(x)\) 的系数,将单位根的倒数 \(w_n^0,w_n^{-1},w_n^{-2},\cdots,w_n^{-n+1}\) 代入以后,得到的每个数再除以 \(n\),就是 \(a(x)\) 的各项系数。

证明

\((b_0,b_1,\cdots,b_{n-1})\)\(A(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^{n-1}\) 离散傅里叶变换的结果。

\(B(x) = b_0+b_1x+\cdots+b_nx^{n-1}\),然后将那几个单位根的倒数代入得到一个新的离散傅里叶变换结果 \((c_0,c_1,\cdots,c_{n-1})\)

\[\begin{aligned}c_k&=\sum_{i=0}^{n-1}y_i(w_n^{-k})^i\\&=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(w_n^i)^j)(w_n^k)^i\\&=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(w_n^{j-k})^i) \end{aligned}\]

发现当且仅当 \(j=k\) 时,后面式子的值为 \(n\),反之,后面的式子值为 \(0\)

那么就有了 \(a_i=\frac{c_i}{n}\)

然后这样就是 \(O(n^2)\) 的了,没啥用啊。

这是因为傅里叶爷爷 (1768年3月21日~1830年5月16日) 没有见过计算机(世界上第一台通用计算机“ENIAC”于1946年2月14日在美国宾夕法尼亚大学诞生),所以他不需要考虑时间复杂度,但是后人就要优化了。

考虑一个多项式 \(A(x)=a_0+a_1x+a_2x^2+\cdots+a^{n-1}x^{n-1}\) 要求离散傅里叶变换,将一个 \(w_n\) 代入。

\(A(x)\) 的每一项按照下标奇偶性分组,设 \(A_1(x)=a_0+a_2x+\cdots+a^{n-2}x^{\frac{n}{2}-1},A_2(x) = a_1 + a_3x + \cdots + a_{n - 1}x^{\frac{n}{2} - 1}\)

显然有 \(A(x) = A_1(x^2)+xA_2(x^2)\)

\(w_n^k(k<\frac{n}{2})\) 代入,有 \(A(w_n^k)=A_1(w_n^{2k})+w_n^kA_2(w_n^{2k})=A_1(w_{\frac{n}{2}}^k)+w_{n}^kA_2(w_{\frac{n}{2}}^k)\)
那么将 \(w_n^{k+\frac{n}{2}}\) 代入,最后得到一个 \(A_1(w_n^{2k})-w_n^kA_2(w_n^{2k})\)

发现这两个式子只有一个常数项不同,那么求第一个式子可以顺便将第二个式子求出来,因为第一个式子取遍 \([0,\frac{n}{2}-1]\),第二个式子取遍 \([\frac{n}{2},n]\),所以将原问题缩小了一半,而且缩小后的问题符合原问题性质,然后就这么递归分治下去即可。

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

【模板】多项式乘法(FFT)

递归实现,跑的好像不是很慢。

个人感觉递归实现便于理解,可以先打递归实现试试水。

code
#include<bits/stdc++.h>
using namespace std;
#define rep(i,s,t,p) for(int i = s;i <= t;i += p)
#define drep(i,s,t,p) for(int i = s;i >= t;i -= p)
#ifdef LOCALauto I = freopen("in.in","r",stdin),O = freopen("out.out","w",stdout);
#elseauto I = stdin,O = stdout;
#endif
using ll = long long;using ull = unsigned long long;
using db = double;using ldb = long double;
using comp = complex<db>;
const db pi = acos(-1);
const int N = 4e6 + 10;
int n,m;comp a[N],b[N];
void fft(comp *a,int n,int inv){if(n == 1) return;//最后了,直接 return 就好啦。int m = n>>1;comp a1[m+5],a2[m+5];rep(i,0,n,2) a1[i>>1] = a[i],a2[i>>1] = a[i+1];//奇偶分组fft(a1,m,inv);fft(a2,m,inv);comp W = comp(cos(2.0*pi/n),inv*sin(2.0*pi/n)),w = comp(1,0);//单位根。for(int i = 0;i < m;++i,w = w*W){a[i] = a1[i] + w*a2[i];//求A1。a[i+m] = a1[i] - w*a2[i];//求A2。}
}
signed main(){cin.tie(nullptr)->sync_with_stdio(false);cin>>n>>m;rep(i,0,n,1) cin>>a[i];rep(i,0,m,1) cin>>b[i];int lim = 1;while(lim <= n+m) lim <<= 1;//凑成 2^nfft(a,lim,1);fft(b,lim,1);//转成点值表示rep(i,0,lim,1) a[i] = a[i]*b[i];//点值乘fft(a,lim,-1);//转成系数表示rep(i,0,n+m,1) cout<<(int)(a[i].real()/lim+0.5)<<' ';
}

但是相比于常写的迭代写法来说,这种写法还是比较慢。(在洛谷这道题,最后一个数据点,递归跑了 1000ms,迭代跑了 600多ms)

那么如何写成迭代法呢?

考虑最后进行操作的序列,发现其实最后操作的序列就是将下标二进制位翻转排序。

手动模拟一下应该很好理解,就是先以二进制位下第零位奇偶分组,然后每组中以第一位奇偶分组,以此类推。

那么就可以知道每个数最后应该在的位置,这个可以递推出来,假如当前位置为 \(i\),那么它会在 (to[i>>1]>>1)|((i&1)<<(ct-1)),其中 \(to_i\) 表示第 \(i\) 个数应该在的位置,显然有 i=to[to[i]]\(ct\) 是二进制位数,即 \(n=2^{ct}\)

然后就可以递推出结果了。

code
#include<bits/stdc++.h>
using namespace std;
#define rep(i,s,t,p) for(int i = s;i <= t;i += p)
#define drep(i,s,t,p) for(int i = s;i >= t;i -= p)
#ifdef LOCALauto I = freopen("in.in","r",stdin),O = freopen("out.out","w",stdout);
#elseauto I = stdin,O = stdout;
#endif
using ll = long long;using ull = unsigned long long;
using db = double;using ldb = long double;
using comp = complex<db>;
const int N = 4e6 + 10;
const db pi = acos(-1);
int n,m,to[N],ct;
comp a[N],b[N];
void fft(comp *a,int n,int type){rep(i,0,n-1,1) if(i < to[i]) swap(a[i],a[to[i]]);for(int mid = 1;mid < n;mid <<= 1){comp W = comp(cos(pi/mid),type*sin(pi/mid));for(int Res = mid<<1,j = 0;j < n;j += Res){comp w = comp(1,0);for(int k = 0;k < mid; ++k,w = w*W){comp x = a[j+k],y = w*a[j+mid+k];a[j+k] = x + y;a[j+mid+k] = x-y;}}}
}
signed main(){cin.tie(nullptr)->sync_with_stdio(false);cin>>n>>m;rep(i,0,n,1) cin>>a[i];rep(i,0,m,1) cin>>b[i];int lim = 1;while(lim <= (n+m)) lim <<= 1,ct++;rep(i,0,lim-1,1) to[i] = (to[i>>1]>>1)|((i&1)<<(ct-1));fft(a,lim,1);fft(b,lim,1);rep(i,0,lim,1) a[i] = a[i]*b[i];fft(a,lim,-1);rep(i,0,n+m,1) cout<<(int)(a[i].real()/lim+0.5)<<' ';
}

实际应用中还可以预处理单位根,但我写丑了,跑的还没不预处理的快(甚至不如递归

code
#include<bits/stdc++.h>
using namespace std;
#define rep(i,s,t,p) for(int i = s;i <= t;i += p)
#define drep(i,s,t,p) for(int i = s;i >= t;i -= p)
#ifdef LOCALauto I = freopen("in.in","r",stdin),O = freopen("out.out","w",stdout);
#elseauto I = stdin,O = stdout;
#endif
using ll = long long;using ull = unsigned long long;
using db = double;using ldb = long double;
using comp = complex<db>;
const int N = 4e6 + 10;
const db pi = acos(-1);
int n,m,to[N],ct;
comp a[N],b[N],urt[N],iurt[N];
void fft(comp *a,comp *urt,int n){rep(i,0,n-1,1) if(i < to[i]) swap(a[i],a[to[i]]);for(int mid = 1;mid < n;mid <<= 1){for(int Res = mid<<1,j = 0;j < n;j += Res){for(int k = 0;k < mid; ++k){comp x = a[j+k],y = urt[n/mid*k]*a[j+mid+k];a[j+k] = x + y;a[j+mid+k] = x-y;}}}
}
signed main(){cin.tie(nullptr)->sync_with_stdio(false);cin>>n>>m;rep(i,0,n,1) cin>>a[i];rep(i,0,m,1) cin>>b[i];int lim = 1;while(lim <= (n+m)) lim <<= 1,ct++;comp W = comp(cos(pi/lim),sin(pi/lim));urt[0] = comp(1,0);rep(i,1,lim,1) urt[i] = urt[i-1]*W;W = conj(W);iurt[0] = comp(1,0);rep(i,1,lim,1) iurt[i] = iurt[i-1]*W;rep(i,0,lim-1,1) to[i] = (to[i>>1]>>1)|((i&1)<<(ct-1));fft(a,urt,lim);fft(b,urt,lim);rep(i,0,lim,1) a[i] = a[i]*b[i];fft(a,iurt,lim);rep(i,0,n+m,1) cout<<(int)(a[i].real()/lim+0.5)<<' ';
}

例题:【模板】高精度乘法 | A*B Problem 升级版

没有压位,和高精一样写即可,就是乘法由暴力乘换成了 fft。

code
#include<bits/stdc++.h>
using namespace std;
#define rep(i,s,t,p) for(int i = s;i <= t;i += p)
#define drep(i,s,t,p) for(int i = s;i >= t;i -= p)
#ifdef LOCALauto I = freopen("in.in","r",stdin),O = freopen("out.out","w",stdout);
#elseauto I = stdin,O = stdout;
#endif
using ll = long long;using ull = unsigned long long;
using db = double;using ldb = long double;
using comp = complex<db>;
const db pi = acos(-1);
const int N = 4e6 + 10;
string s1,s2;
int n,m,ct,to[N],ans[N];comp a[N],b[N];
void fft(comp *a,int type,int n){rep(i,0,n,1) if(i < to[i]) swap(a[i],a[to[i]]);for(int mid = 1;mid < n;mid <<= 1){comp W = comp(cos(pi/mid),type*sin(pi/mid));for(int Res = mid<<1,j = 0;j < n;j += Res){comp w = comp(1,0);for(int k = 0;k < mid; ++k,w = w*W){auto x = a[j+k],y = w*a[j+mid+k];a[j+k] = x+y;a[j+mid+k] = x-y;}}}
}
signed main(){cin.tie(nullptr)->sync_with_stdio(false);cin>>s1>>s2;n = s1.size() - 1,m = s2.size() - 1;reverse(s1.begin(),s1.end());reverse(s2.begin(),s2.end());rep(i,0,n,1) a[i] = comp(s1[i]-'0',0);rep(i,0,m,1) b[i] = comp(s2[i]-'0',0);int lim = 1;while(lim <= (n+m)) lim <<= 1,ct++;rep(i,0,lim,1) to[i] = (to[i>>1]>>1)|((i&1)<<(ct-1));fft(a,1,lim);fft(b,1,lim);rep(i,0,lim,1) a[i] = a[i]*b[i];fft(a,-1,lim);rep(i,0,lim,1) ans[i] = (int)(a[i].real()/lim+0.5);rep(i,0,lim,1) ans[i+1] += ans[i]/10,ans[i] %= 10;int now = lim;while(!ans[now]) now--;drep(i,now,0,1) cout<<ans[i];cout<<'\n';
}

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.hqwc.cn/news/878473.html

如若内容造成侵权/违法违规/事实不符,请联系编程知识网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

GPUStack:一个开源的GPU集群管理和自动化部署大模型的LLM服务平台

正如我们以前文章中说过,在过去近两年中,AI发展速度超过任何历史时期,基于大模型的各类应用已经渗透到了各行各业中,很多企业都在积极探索如何利用大模型提高公司运营管理的能效。 阿里云 CTO 周靖人也说过““当下企业应用大模型存在三种范式:一是对大模型开箱即用,二是…

Deskflow:一个能在多个设备之间共享鼠标键盘的工具

有多台设备的同学肯定都感受过切换鼠标和键盘的痛苦。比如家里有个台式机配了键盘和鼠标,工作有台笔记本,但是如果想用家里的键盘和鼠标操作笔记本那么除了插拔线(蓝牙也一样),好像没有其他特别的办法。 以为我个人为例,插拔切换鼠标键盘,最头疼的就是以下两个点: 线缆…

DeepSeek + Dify :零成本搭建企业级本地私有化知识库保姆级喂饭教程

最近,DeepSeek大火,想必大家都有所耳闻,各路媒体从各个方面报道了DeepSeek这家神秘的公司的各方面消息,这家低调的技术公司用一组硬核数据回应了所有关注: 千亿参数规模下实现0.5元/百万tokens的API调用成本,91.5%的中文基准测试得分,推理效率较传统架构提升5倍。 DeepS…

[Tools] Vite intro

Overview为什么选Vite: https://cn.vite.dev/guide/why.htmlesbuild, Rollup: https://cn.vite.dev/guide/why.html#why-bundle-for-production Quick start 1. Start a new project pnpm init 2. installl Vite, Typescript pnpm add vite typescript -D 3. create index.ht…

在做同城多活方案中如何实现机房之间的数据同步?

一、前言 现在是数据的时代,如何发挥数据的价值,让系统具备更多的数据处理能力。如何完善响应的数据架构。多机房建设是其中的解决策略。 二、背景 在业务初期,考虑到投入成本,很多公司通常只用一个机房提供服务。但随着业务发展,流量不断增加,我们对服务的响应速度和可用…

在 Ubuntu 22.04 上运行 Filebeat 7.10.2

环境 操作系统:阿里云 Ubuntu 22.04.3 LTS (GNU/Linux 5.15.0-83-generic x86_64) 软件版本:Filebeat 7.10.2 用户:root 运行下载从这里下载 filebeat 7.10.2。配置简单配置一下 filebeat.yml,从标准输入采集,写入到标准输出 : filebeat.inputs: - type: stdinoutput.con…

Dify AWS:0代码搭建「AI翻译中台」

0 前言 基于Dify现有能力,已能对不少业务场景提供帮助,但对一些特定诉求,还要借助其扩展机制,本文利用翻译场景举例详细说明。 1 翻译场景复杂性分析 翻译是从简单到复杂各级都存在的场景,比较简单的翻译可能一句简单 Prompt,但对复杂、效果要求较高翻译场景,可能需要一…

25.2.3小记

抽象(abstract) 1.抽象函数不需要大括号 2.抽象的函数只能包括在抽象的类中(或者说只要一个类里有一个抽象函数,那整个类就是抽象的,因为若非抽象则定义对应类的变量之后本应该可以调用name.work类似的函数,但抽象函数因为没有body所以无法调用)其中定义的变量可以管理任…

0x80070035错误怎么解决?Win11/Win10访问NAS smba共享文件夹提示无法找到路径[Path not found]

1. 问题 Win11/Win10访问NAS smba共享文件夹提示无法找到路径[Path not found]2. 解决办法 2.1 设置网络专用网络2.2 设置高级共享2.3 设置组策略-启用-不安全的来宾登录2.4 设置组策略-Microsoft网络客户端:对通信进行数字签名(始终)-改为禁用3. 验证 现在再去访问,你会发现…

DeepSeek,你是懂.NET的!

这两天火爆出圈的话题,除了过年,那一定是DeepSeek!你是否也被刷屏了?DeepSeek 是什么DeepSeek是一款由国内人工智能公司研发的大型语言模型,拥有强大的自然语言处理能力,能够理解并回答问题,还能辅助写代码、整理资料和解决复杂的数学问题。与OpenAI开发的ChatGPT相比,De…