浅谈算法——多项式乘法相关
从多项式乘法到FFT
这一段大部分是复制以前我写的这篇博客：
https://blog.csdn.net/wang3312362136/article/details/79510933
这篇博客有详细的代码，但是阅读体验布星，因此我把它的大部分证明过程复制下来。
如果会fft，请跳过这一段。
多项式乘法
我们知道，多项式可以表示成：
A=\sum_{i=0}^{n}a_ix^i的形式。
对于两个多项式A(x)A(x)A(x)和B(x)B(x)B(x)，我们可以计算乘积A⋅BA⋅BA\cdot B：
A\cdot B=\sum_{i=0}^{|A|}\sum_{j=0}^{|B|}a_ib_jx^{i+j}
但是，这样算是O(|A|×|B|)O(|A|×|B|)O(|A|\times |B|)的，太慢了，怎么办？
我们需要换一条思路。
首先，我们得知道一个东西：多项式的点值表示法。
我们把上面的称为多项式的系数表示法，而点值表示法就是：
若AAA多项式的次数为n" role="presentation">nnn，则任取nnn个不相同的x0,x1,⋯,xn" role="presentation">x0,x1,⋯,xnx0,x1,⋯,xnx_0,x_1,\cdots,x_n，求出AAA多项式的A(x0),A(x1),⋯,A(xn)" role="presentation">A(x0),A(x1),⋯,A(xn)A(x0),A(x1),⋯,A(xn)A(x_0),A(x_1),\cdots,A(x_n)。记为：
显然，一个有n+1n+1n+1个点的点对唯一表示一个nnn次多项式。
对于一个点值表示法下多项式
和
它们的乘积是
可以看出点值表示法的多项式相乘是O(n)O(n)O(n)的。
等等，我们好像找到了一个突破口！
为啥不把原来的系数表示法下的多项式转化成点值表示法呢？
仔细想一想：系数表示法与点值表示法互相转换，这个步骤好像是O(n2)O(n2)O(n^2)的。
而FFT（快速傅里叶变换）就是为了优化这个O(n2)O(n2)O(n^2)。
PS：对于O(n2)O(n2)O(n^2)的点值表示法转化成系数表示法可以看百度百科中关于插值法的介绍。
FFT（快速傅里叶变换）
如果未特别说明，那么下面的多项式次数将是2k−12k−12^k-1的形式。
如果不是关键部分的公式或定理，不提供证明，自己出门右转百度。
首先介绍两个概念：
DFT（离散傅里叶变换）是将多项式由系数表示法转化成点值表示法；
IDFT（离散傅里叶逆变换）是将多项式由点值表示法转化成系数表示法；
而FFT就是上述两种变换的优化。
DFT部分
前置技能：
下面的内容将会提到复数，不会的可以参考百度百科中关于复数的介绍；
为了介绍FFT中的DFT部分，首先要介绍的是一个概念：单位根。
单位根：若有
z^n=1此时将 zzz称为n" role="presentation">nnn次单位根。
若有z∈ℝz∈Rz\in \mathbb R，显然，zzz可以等于1" role="presentation">111，如果nnn是偶数，则z" role="presentation">zzz还可以等于−1−1-1。
我们把范围扩大到z∈ℂz∈Cz \in \mathbb C，那么，我们可以得到nnn个复数，它们将复平面上的单位圆等分成n" role="presentation">nnn份。
为了表示nnn次单位根，我们引入一个公式。
欧拉公式：
e^{in}=\cos n+i\sin n
如果我们令：
\omega_n=e^{2\pi i/n}
那么，nnn次单位根就可以表示成ωn0,ωn1,⋯,ωnn−1" role="presentation">ω0n,ω1n,⋯,ωn−1nωn0,ωn1,⋯,ωnn−1\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}，它们的nnn次方显然都是1" role="presentation">111。
下面是关于ωnωn\omega_n的两条性质：（都是在nnn为偶数的情况下）
\omega_n^{n/2}=e^{(2\pi i/n)\cdot (n/2)}=e^{\pi i}=cos \pi+i\sin \pi=-1\tag{1}
\omega_n^2=e^{2\cdot 2\pi i/n}=\omega_{n/2}\tag{2}
下面，我们进入正题：DFT的求法。
在这里，我们令多项式次数为n−1n−1n-1，那么我们可以用点值表示成
额……这时间复杂度好像并没有减少……
别急，我们来看A(ωkn)A(ωnk)A(\omega_n^k)能够表示成什么。
\begin{align}A(\omega_n^k)=&\sum_{i=0}^{n-1}a_i\omega_n^{ki} \nonumber\tag{3}\\=&\sum_{i=0}^{n/2-1}a_{2i}\omega_n^{2ki}+\sum^{n/2-1}_{i=0}a_{2i+1}\omega_n^{2ki+k}\nonumber\tag{4}\\=&\sum_{i=0}^{n/2-1}a_{2i}\omega_{n/2}^{ki}+\sum_{i=0}^{n/2-1}a_{2i+1}\omega_{n/2}^{ki}\omega_{n}^k\nonumber\tag{5}\\=&\sum_{i=0}^{n/2-1}a_{2i}\omega_{n/2}^{ki}+\omega_n\sum_{i=0}^{n/2-1}a_{2i+1}\omega_{n/2}^{ki}\nonumber\tag{6}\end{align}
我们来分别分析以下这几个神奇的步骤。
(3)(3)(3)步骤就是将ωknωnk\omega_n^k带入原来的AAA多项式。
(4)" role="presentation">(4)(4)(4)步骤就是将原多项式拆成两个部分，按奇偶分类。
(5)(5)(5)步骤用到了上面提到的性质(2)(2)(2)。
(6)(6)(6)步骤就是上面式子的后半部分提出公因数。
有了这个等式，我们就可以分治+递归解决DFT了。
算法步骤：
- 对当前的多项式（一个数组）系数进行奇偶分类；
- 递归算出偶数部分的数组的anseanseans_e和奇数部分的数组的ansoansoans_o；
- 这个多项式的ans=anse+ωnansoans=anse+ωnansoans=ans_e+\omega_nans_o
但是这个的常数好像很大啊？能不能减少一点呢？
上面的性质(1)(1)(1)给了我们提示：
\omega_n^{n/2+k}=\omega_n^{n/2}\cdot \omega_n^{k}=-\omega_n^k
在算k<n2k<n2k时，可以顺便把k≥n2k≥n2k\geq\frac{n}{2}的情况也算出来。
常数减小了一半！但是还是很大啊！
递归版的程序一般比非递归版慢，为啥不用非递归版呢？
算法核心就是奇偶分类，分来分去最后分到了哪里？我们来研究研究。
显然，一个序列原来是0,1,2,3,4,5,6,70,1,2,3,4,5,6,70,1,2,3,4,5,6,7，最终变成0,4,2,6,1,5,3,70,4,2,6,1,5,3,70,4,2,6,1,5,3,7。
把它们的二进制列出来：
000,001,010,011,100,101,110,111\\000,100,010,110,001,101,011,111
其中，上面是位置，下面是这个位置对应的数。
把上面的数翻转，好像就是下面的数！
没错，只需要计算一下每个数的二进制翻转后的结果，就能得到一个数最终对应的位置，也就能实现非递归版了。
IDFT部分
回顾上面的DFT部分，仔细思考一下，它本质就是在求：
\begin{cases}a_0(\omega_n^0)^0+a_1(\omega_n^0)^1+\cdots+a_{n-1}(\omega_n^0)^{n-1}=A(\omega_n^0)\\a_0(\omega_n^1)^0+a_1(\omega_n^1)^1+\cdots+a_{n-1}(\omega_n^1)^{n-1}=A(\omega_n^1)\\\cdots\\a_0(\omega_n^{n-1})^0+a_1(\omega_n^{n-1})^1+\cdots+a_{n-1}(\omega_n^{n-1})^{n-1}=A(\omega_n^{n-1})\end{cases}
其中，给定了a0,a1,⋯,an−1a0,a1,⋯,an−1a_0,a_1,\cdots, a_{n-1}以及ω0n,ω1n,⋯,ωn−1nωn0,ωn1,⋯,ωnn−1\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}，求A(ω0n),A(ω1n),⋯,A(ωn−1n)A(ωn0),A(ωn1),⋯,A(ωnn−1)A(\omega_n^0),A(\omega_n^1),\cdots,A(\omega_n^{n-1})的值。
用矩阵表示如下：
\begin{bmatrix}(\omega_n^0)^0&(\omega_n^0)^1&\cdots&(\omega_n^{0})^{n-1}\\(\omega_n^1)^0&(\omega_n^1)^1&\cdots&(\omega_n^1)^{n-1}\\\vdots&\vdots&\ddots&\vdots\\(\omega_n^{n-1})^0&(\omega_n^{n-1})^1&\cdots&(\omega_n^{n-1})^{n-1}\end{bmatrix}\begin{bmatrix}a_0\\a_1\\\vdots\\a_{n-1}\end{bmatrix}=\begin{bmatrix}A(\omega_n^0)\\A(\omega_n^1)\\\vdots\\A(\omega_n^{n-1})\end{bmatrix}\tag{7}
我们令：
V=\begin{bmatrix}(\omega_n^0)^0&(\omega_n^0)^1&\cdots&(\omega_n^{0})^{n-1}\\(\omega_n^1)^0&(\omega_n^1)^1&\cdots&(\omega_n^1)^{n-1}\\\vdots&\vdots&\ddots&\vdots\\(\omega_n^{n-1})^0&(\omega_n^{n-1})^1&\cdots&(\omega_n^{n-1})^{n-1}\end{bmatrix}
那么IDFT的本质就是求VVV矩阵的逆矩阵。
普通的矩阵求逆显然无法快速求出，这必须依靠人类智慧大力猜想逆矩阵。
人类智慧还真求出来了，考虑下面这个矩阵：
D=\begin{bmatrix}(\omega_n^{-0})^0&(\omega_n^{-0})^1&\cdots&(\omega_n^{-0})^{n-1}\\(\omega_n^{-1})^0&(\omega_n^{-1})^1&\cdots&(\omega_n^{-1})^{n-1}\\\vdots&\vdots&\ddots&\vdots\\(\omega_n^{-(n-1)})^0&(\omega_n^{-(n-1)})^1&\cdots&(\omega_n^{-(n-1)})^{n-1}\end{bmatrix}
那么我们令E=D⋅VE=D⋅VE=D\cdot V，则：
\begin{align}E_{i,j}=&\sum_{k=0}^{n-1}D_{i,k}V_{k,j}\nonumber\\=&\sum_{k=0}^{n-1}(\omega_n^{-i})^k(\omega_n^k)^j\nonumber\\=&\sum_{k=0}^{n-1}\omega_n^{k(j-i)}\nonumber\end{align}
显然：
E_{i,j}=\begin{cases}0(i\not= j)\\n(i=j)\end{cases}
因此
\frac{1}{n}D\cdot V=\frac{1}{n}E=I_n
(7)(7)(7)式两边同时左乘1nD1nD\frac{1}{n}D，可得
\begin{bmatrix}a_0\\a_1\\\vdots\\a_{n-1}\end{bmatrix}=\frac{1}{n}\begin{bmatrix}(\omega_n^{-0})^0&(\omega_n^{-0})^1&\cdots&(\omega_n^{-0})^{n-1}\\(\omega_n^{-1})^0&(\omega_n^{-1})^1&\cdots&(\omega_n^{-1})^{n-1}\\\vdots&\vdots&\ddots&\vdots\\(\omega_n^{-(n-1)})^0&(\omega_n^{-(n-1)})^1&\cdots&(\omega_n^{-(n-1)})^{n-1}\end{bmatrix}\begin{bmatrix}A(\omega_n^0)\\A(\omega_n^1)\\\vdots\\A(\omega_n^{n-1})\end{bmatrix}
这就相当于把DFT中ωknωnk\omega_n^k都换成ω−knωn−k\omega_{n}^{-k}。
IDFT过程中的小技巧
如果IDFT中逆矩阵使用的是VVV（使用的是DFT的过程），要得到IDFT的正确结果只需要执行一遍std::reverse(a+1,a+n)，这是很容易证明的。
NTT（快速数论变换）
FFT要求采用复数和实数运算，如果数很大很容易产生精度误差。如果是模意义下呢？
前置技能
FFT，初级数论内容，原根。
一点微小的变化？
考虑原根的性质。
性质(1)" role="presentation">(1)(1)(1) 循环：gn+xn=gxn(modn)gnn+x=gnx(modn)g_n^{n+x}=g_n^x \pmod{n}
性质(2)(2)(2) 折半：gxn=gx/2n/2(modm)gnx=gn/2x/2(modm)g_n^x=g_{n/2}^{x/2}\pmod{m}
性质(3)(3)(3) ∑n−1i=0gin=0(modm)∑i=0n−1gni=0(modm)\sum_{i=0}^{n-1} g_n^i=0\pmod{m}
可以发现原根的性质和单位根的性质差不多，因此可以用来做模意义下的FFT，就是NTT。
或者说，NTT就是把FFT的ωω\omega替换成了ggg。
由于使用的是gim−1" role="presentation">gm−1igim−1g_i^{m-1}，其中mmm为模数，i" role="presentation">iii为2k2k2^k形式，因此模数是x⋅2k+1x⋅2k+1x\cdot 2^k+1的形式（NTT模数）。
常见的NTT模数：
998244353，原根是3
1004535809，原根也是3
FWT（快速沃尔什变换）
上面求的都是多项式卷积
F(n)=\sum_{i=0}^n a(i)b(n-i)
或者换种方式
F(n)=\sum_{i,j,i+j=n}a(i)b(j)
但是如果是下面的方式
F(n)=\sum_{i,j,i\circ j=n} a(i)b(j)
其中 ∘∘\circ是种位运算呢？
依然考虑转化成点值表达式再相乘。类比FFT，考虑用一种变换tftftf和逆变换utfutfutf，只要tf(A×B)=tf(A)×tf(B)tf(A×B)=tf(A)×tf(B)tf(A\times B)=tf(A)\times tf(B)，那么就可以得到A×B=utf(tf(A)×tf(B))A×B=utf(tf(A)×tf(B))A\times B=utf(tf(A)\times tf(B))
下面我们分类讨论∘∘\circ的位运算种类及tf,utftf,utftf,utf的种类。
集合交卷积（and）
集合并卷积（or）
集合对称差卷积（xor）
子集卷积
其他的方式
多项式求逆
前置技能
多项式乘法（FFT/NTT）
多项式的度
一个多项式的度就是这个多项式最高此项的次数，例如4x2+3x+14x2+3x+14x^2+3x+1的度为222。
多项式的逆元
现在有两个多项式F,G" role="presentation">F,GF,GF,G，如果
F\times G=1 \pmod{x^n}\tag{1}
并且 GGG的度小于等于F" role="presentation">FFF的度，那么 GGG称为F" role="presentation">FFF在模 xnxnx^n意义下的逆元。
怎么求？假设已经求出了FFF在模x⌈n/2⌉" role="presentation">x⌈n/2⌉x⌈n/2⌉x^{\lceil n/2\rceil}下的逆元，为G′G′G'。
由定义得到
F\times G'=1 \pmod{x^{\lceil n/2\rceil}}\tag{2}
将(1)(1)(1)放到模x⌈n/2⌉x⌈n/2⌉x^{\lceil n/2\rceil}意义下
F\times G=1\pmod{x^{\lceil n/2\rceil}}\tag{3}
(2)−(3)(2)−(3)(2)-(3)得到
F\times G'-F\times G=0 \pmod{x^{\lceil n/2\rceil}}
两边平方，得到
F^2\times G^2-2\times F^2\times G'\times G+F^2\times G'^2=0\pmod{x^n}
同时除以 FFF并处理modxn" role="presentation">modxnmodxn\bmod x^{n}意义下为 111的项
G=2G'-FG'^2\pmod {x^n}
就是说，只要求出了模 x⌈n/2⌉x⌈n/2⌉x^{\lceil n/2\rceil}意义下的逆元，就可以求出模 xnxnx^n意义下的逆元。
浅谈算法——多项式乘法相关相关推荐
- 浅谈用户密码保护与相关技术
浅谈用户密码保护与相关技术(上) 一. 全文涉及 上篇:哈希,彩虹表 下篇:加盐加密,慢哈希,非对称加密与HTTPS 二. 主题引入 2011年12月21日,CSDN后台数据库被黑客恶意发布到互联 ...
- 浅谈算法和数据结构: 五 优先级队列与堆排序
原文:浅谈算法和数据结构: 五 优先级队列与堆排序 在很多应用中,我们通常需要按照优先级情况对待处理对象进行处理,比如首先处理优先级最高的对象,然后处理次高的对象.最简单的一个例子就是,在手机上玩游戏 ...
- 算法工程师_浅谈算法工程师的职业定位与发展
随着大数据和以深度学习为代表的人工智能技术的飞速发展,算法工程师这个职业逐渐成为国内互联网行业的标配.2016年3月,谷歌旗下DeepMind公司的围棋程序"AlphaGo"战胜职 ...
- 【学习笔记】浅谈广义矩阵乘法——动态DP
文章目录 广义矩阵乘法 动态DP 例题:洛谷4719 以下内容是本人做题经验,如有雷同,纯属抄袭:如有不对,纯属不懂,还请指正 广义矩阵乘法 众所周知,矩阵满足乘法交换律,前一个矩阵的列必须是后一个矩 ...
- 从石头剪刀布浅谈算法的作用
刚开始学习C语言的时候,常常听到前辈说,C语言的核心就是算法.但是对于小白来说,常常一脸懵逼,搞不懂啥叫算法?算法有什么用?我的if-else语句照样可以走天下.但是作为小白来说虽然不懂但是也不敢问, ...
- 浅谈Java中类的相关内容
目录 阅前说明 面向对象思想 面向对象思想概述 类和对象 代码示例 对象内存图 继承 概述 定义 优点 格式 代码示例 super和this super和this的含义 用法 抽象类 概述 抽象类的使 ...
- 浅谈算法和数据结构: 七 二叉查找树
前文介绍了符号表的两种实现,无序链表和有序数组,无序链表在插入的时候具有较高的灵活性,而有序数组在查找时具有较高的效率,本文介绍的二叉查找树(Binary Search Tree,BST)这一数据结构 ...
- 浅谈算法（简单算法）
前言 算法(Algorithm)是指解题方案的准确而完整的描述,是一系列解决问题的清晰指令,算法代表着用系统的方法描述解决问题的策略机制.也就是说,能够对一定规范的输入,在有限时间内获得所要求的输出. ...
- 浅谈Android之SurfaceFlinger相关介绍（一）
SurfaceFlinger是GUI的核心,以系统服务的形式存在,负责将所有App的图形数据按照Z Order顺序混合并输出到FrameBuffer. 根据图中描述,从下到上依次介绍: 1) 这里的 ...
最新文章
- windows7升级安装之初体验
- linux 进程内存开销,linux下查看最消耗CPU、内存的进程
- 【Python】pyCryptodome模块实现AES加密、解密
- 解读中国互联网：局部领先、快进的数字化发展
- Eclipse安装后启动出现error:could not create the java machine.
- latex 论文绘图： 图像文字重叠
- 用python简单代码做一个计算器
- ceph-cache-tier
- 这7个摄影构图技巧，可能会帮你拍出好看照片！你学会了吗？
- Java 生成数字证书系列（四）生成数字证书（续）
- 第52届格莱美大奖完全获奖名单
- Could not find or load main class org.apache.hadoop.mapreduce.v2.app.MRAppMaster
- html返回首页页面代码,后台返回的HTML整个页面代码打开方法
- 回顾小米公司的成功过程，用五个层次的问题阐述，小米的成功基础、小米的爆品特点、小米生态链模式的根本原因。
- yarn create umi 报错问题
- R语言读写中文编码方式
- 自适应波束形成（四）——Frost波束形成1
- 笔记本120赫兹输出html,120Hz显示器vs.60Hz显示器盲测
- 学专业计算机可以当游戏主播吗,一个专业的游戏主播需要什么配置的电脑
- SEXTANTE中调用任意C++控制台程序的简单例子