FFT求解卷积

最新推荐文章于 2023-06-06 23:41:22 发布

置顶兔蠢蠢

最新推荐文章于 2023-06-06 23:41:22 发布

阅读量1.4k

点赞数 1

分类专栏：数论数学 FFT

本文链接：https://blog.csdn.net/ftx456789/article/details/79172848

版权

数论同时被 3 个专栏收录

33 篇文章 0 订阅

订阅专栏

数学

28 篇文章 0 订阅

订阅专栏

FFT

4 篇文章 0 订阅

订阅专栏

卷积

已知两个多项式
$A(x)=a_{0}+a_{1}x^1+a_{2}x^2+a_{3}x^3+\cdots +a_{n}x^n$
$B(x)=b_{0}+b_{1}x^1+b_{2}x^2+b_{3}x^3+\cdots +b_{m}x^m$
$C (x) = A (x) B (x)$
那么多项式C(x)得到的是一个n+m-1的多项式

$C(x)=c_{0}+c_{1}x^1+c_{2}x^2++c_{3}x^3+\cdots +c_{n+m-1}x^{n+m-1}$

$C(x)=\sum_{i=0}^{n+m-1}c_{i}x^{i}$

其中 $c_{i}=\sum_{j=0}^{i}a_{j}b_{i-j}$

我们可以把多项式A(x),B(x),C(x)的系数以一个行向量来表示

$\vec{a}=\left ( a_{0},a_{1},a_{2},\cdots ,a_{n} \right )$

$\vec{b}=\left ( b_{0},b_{1},b_{2},\cdots ,b_{m} \right )$

$\vec{c}=\left ( c_{0},c_{1},c_{2},\cdots ,c_{n+m-1} \right )$

通过向量 $\vec{a}$ 和 $\vec{b}$ 求解出 $\vec{c}$ 的操作我们称为卷积

$\vec{c}=\vec{a}\bigotimes \vec{b}$ ,$\bigotimes $为卷积符号

相关算法

求卷积的算法直接模拟的话时间复杂度是 $O(n^{2})$ 的,这个时候我们有一个算法来优化实现 $O (n l o g (n))$ 的时间复杂度

就像我们前面讲的多项式的系数可以表示成行向量的形式，一个多项式也可以用这个多项式上的一些点值来表示如

比如多项式 $y = A (x)$ 可以用 $\left \{ (x_{0},y_{0}),(x_{1},y_{1}),(x_{2},y_{2}),\cdots ,(x_{n},y_{n}) \right \}$ 这样的n个点来表示

举个具体的例子
$A(x)=x^2-2x+3$
$\left \{ (0,3),(1,2),(2,3),(3,6),(4,11) \right \}$

$B(x)=x^2+x+1$
$\left \{ (0,1),(1,3),(2,7),(3,13),(4,21) \right \}$

$C(x)=A(x)B(x)=x^4-x^3+2x^2+3$
$\left \{ (0,3),(1,6),(2,21),(3,78),(4,231) \right \}$

从点值表达式的例子上可以看出,两个多项式的点值表达式的乘积的值与这两个多项式乘完之后的多项式的点值表达式是一致的,那么就有一个想法我们可不可以把多项式转换为点值表达式,然后通过点值表达式的乘积获得原两个多项式乘积后的多项式的点值表达式,再又这个点值表达式来求得多项式呢。这显然是可以的。
把多项式转化为点值表达式的算法为DFT(离散傅里叶变换)
把点值表达式转换为多项式的算法为IDFT(逆DFT)
利用DFT和IDFT来求解多项式的方法叫做FFT
DFT原理
我们前面说了要构建一个点值表达式,点值表达式中要求每一个 $x_{i}$ 都要不同,这个时候我们就可以用到复数的性质了简化运算了
复数我们可以用指数形式来表达
$z = a + i b$
$z=re^{iθ}$
这里r为复数z的模,θ为幅角度数
再回到 $x_{i}$ 的取值,我们可以取 $w^{n}=1$ ,的解作为 $x_{i}$ 的取值,其中 $w^{n}$ 恰是1的n次单位复根,其个数也为n个
即用指数形式表达的话为 $e^{\frac{2\pi ki}{n}}$ , $k=0,1,2,3,4,\cdots ,n-1$ ,即有n个

n次单位复根有一些性质

相消引理
对于任意 $k > 0, d > 0$ ,有 $w_{dn}^{dk}=w_{n}^{k}$

证明: $w_{n}^{k}=e^{\frac{2\pi ki}{n}}$ $w_{dn}^{dk}=e^{\frac{d2\pi ki}{nd}}=e^{\frac{2\pi ki}{n}}$ 即 $w_{dn}^{dk}=w_{n}^{k}$
2. 折半引理
如果有n>0,且n为偶数有 $(w_{n}^{k+\frac{n}{2}})^{2}=(w_{n}^{k})^{2}$

证明: $(e^{\frac{2\pi (k+\frac{n}{2})i}{n}})^{2}$ = $(e^{\frac{2\pi ki}{n}}e^{\frac{2\pi ni}{2n}})^2=(e^{\frac{2\pi ki}{n}})^2$
3. 求和引理
对于 $n > = 1$ 且 $n$ 不能整除 $k$ ,有 $\sum_{i=0}^{n-1}(w_{n}^{k})^i=0$

证明: $\sum_{i=0}^{n-1}(w_{n}^{k})^i$ 是一个等比数列求和所以用等比数列求和公式有

$\sum_{i=0}^{n-1}(w_{n}^{k})^i=\frac{1-(w_{n}^{k})^{n}}{1-w_{n}^{k}}=\frac{1-(w_{n}^{n})^{k}}{1-w_{n}^{k}}=\frac{1-(e^{\frac{2\pi ki}{n}n})^{k}}{1-w_{n}^{k}}=\frac{1-1}{1-w_{n}^{k}}=0$

然后让我们回到 $A (x)$ 来

原来的 $A(x)=\sum_{i=0}^{n}a_{i}x^{i}$ 那么我们现在让A(x)多增加一些项,让他的项数增加到 $len=2^{t}$ , $l e n$ 即为大于n的2的次幂,多余的项我们可以把系数 $a_{i}=0$
为了方便我们令新的 $n = l e n$

我们把 $w_{n}^{k}$ 代入 $A (x)$ ,再将 $A (x)$ 的系数分为奇偶两个式子有

$A_{0}(w_{n}^{k})=a_{0}+a_{2}w_{n}^{k}+a_{4}(w_{n}^{k})^{2}+\cdots +a_{n-2}(w_{n}^{k})^{\frac{n}{2}-1}$

$A_{1}(w_{n}^{k})=a_{1}+a_{3}w_{n}^{k}+a_{5}(w_{n}^{k})^{2}+\cdots +a_{n-1}(w_{n}^{k})^{\frac{n}{2}-1}$

那么原式就有

$A(w_{n}^{k})=A_{0}((w_{n}^{k})^2)+w_{n}^{k}A_{1}((w_{n}^{k})^{2})$

$=A_{0}((e^{\frac{2\pi ki}{n}})^2)+w_{n}^{k}A_{1}((e^{\frac{2\pi ki}{n}})^{2})$

$=A_{0}(e^{\frac{2\pi ki}{\frac{n}{2}}})+w_{n}^{k}A_{1}(e^{\frac{2\pi ki}{\frac{n}{2}}})$

$=A_{0}(w_{\frac{n}{2}}^{k})+w_{n}^{k}A_{1}(w_{\frac{n}{2}}^{k})$

我们再利用折半原理

$A(w_{n}^{k+\frac{n}{2}})=A_{0}((w_{n}^{k+\frac{n}{2}})^2)+w_{n}^{k+\frac{n}{2}}A_{1}((w_{n}^{k+\frac{n}{2}})^{2})$

$=A(w_{n}^{k+\frac{n}{2}})=A_{0}((w_{n}^{k})^2)+w_{n}^{k+\frac{n}{2}}A_{1}((w_{n}^{k})^{2})$

$=A_{0}(w_{\frac{n}{2}}^{k})+w_{n}^{k+\frac{n}{2}}A_{1}(w_{\frac{n}{2}}^{k})$

又有 $w_{n}^{k+\frac{n}{2}}=e^{\frac{2\pi (k+\frac{n}{2})i}{n}}=e^{\frac{2\pi ki}{n}}e^{\pi i}=-e^{\frac{2\pi ki}{n}}=-w_{n}^{k}$

所以原式有

$=A_{0}(w_{\frac{n}{2}}^{k})-w_{n}^{k}A_{1}(w_{\frac{n}{2}}^{k})$

由于 $A(w_{n}^{k+\frac{n}{2}})$ 和 $A(w_{n}^{k})$ 由两个一样的多项式可以求得,这就代表说所代入的单位复根的数目也可以折半了
同样的 $A_{0}(x)$ 和 $A_{1}(x)$ 也可以通过这样的操作再分成奇偶项,这样从迭代下去,这样优化求点值表达式的话时间复杂度就变成 $O (n l o g n)$
具体求法还涉及了蝴蝶操作
求完两个多项式的点值表达式以后我们再把他们相乘,就得到了新的多项式的点值表达式了,那么如何把一个点值表达式再转化为多项式呢
IDFT原理
由点值表达式求解系数可以看作是一个解n元线性方程组的过程,那么我们就得到了一个方程组
$\left\{\begin{matrix} a_{0}(w_{n}^{0})^{0}+a_{1}(w_{n}^{0})^{1}+\cdots +a_{n-1}(w_{n}^{0})^{n-1}=A(w_{n}^{0})\\ a_{0}(w_{n}^{1})^{0}+a_{1}(w_{n}^{1})^{1}+\cdots +a_{n-1}(w_{n}^{1})^{n-1}=A(w_{n}^{1})\\ \vdots\\ a_{0}(w_{n}^{n-1})^{0}+a_{1}(w_{n}^{n-1})^{1}+\cdots +a_{n-1}(w_{n}^{n-1})^{n-1}=A(w_{n}^{n-1}) \end{matrix}\right.$
我们写出方程组的矩阵形式
$\begin{bmatrix} (w_{n}^{0})^{0} & (w_{n}^{0})^{1} &\cdots & (w_{n}^{0})^{n-1}\\ (w_{n}^{1})^{0} & (w_{n}^{1})^{1} &\cdots & (w_{n}^{1})^{n-1}\\ \vdots & \vdots & \ddots & \vdots \\ (w_{n}^{n-1})^{0} & (w_{n}^{n-1})^{1} &\cdots & (w_{n}^{n-1})^{n-1} \end{bmatrix}\begin{bmatrix} a_{0}\\ a_{1}\\ \vdots \\ a_{n-1} \end{bmatrix}=\begin{bmatrix} A(w_{n}^{0})\\ A(w_{n}^{1})\\ \vdots \\ A(w_{n}^{n-1}) \end{bmatrix}$
我们可以求解
$D=\begin{bmatrix} (w_{n}^{0})^{0} & (w_{n}^{0})^{1} &\cdots & (w_{n}^{0})^{n-1}\\ (w_{n}^{1})^{0} & (w_{n}^{1})^{1} &\cdots & (w_{n}^{1})^{n-1}\\ \vdots & \vdots & \ddots & \vdots \\ (w_{n}^{n-1})^{0} & (w_{n}^{n-1})^{1} &\cdots & (w_{n}^{n-1})^{n-1} \end{bmatrix}$ 的逆矩阵
可得
$D^{-1}=\frac{1}{n}\begin{bmatrix} (w_{n}^{-0})^{0} & (w_{n}^{-0})^{1} &\cdots & (w_{n}^{-0})^{n-1}\\ (w_{n}^{-1})^{0} & (w_{n}^{-1})^{1} &\cdots & (w_{n}^{-1})^{n-1}\\ \vdots & \vdots & \ddots & \vdots \\ (w_{n}^{-(n-1)})^{0} & (w_{n}^{-(n-1)})^{1} &\cdots & (w_{n}^{-(n-1)})^{n-1} \end{bmatrix}$
我们把系数矩阵两边同时乘上一个 $D^{-1}$ 有
$\frac{1}{n}\begin{bmatrix} (w_{n}^{-0})^{0} & (w_{n}^{-0})^{1} &\cdots & (w_{n}^{-0})^{n-1}\\ (w_{n}^{-1})^{0} & (w_{n}^{-1})^{1} &\cdots & (w_{n}^{-1})^{n-1}\\ \vdots & \vdots & \ddots & \vdots \\ (w_{n}^{-(n-1)})^{0} & (w_{n}^{-(n-1)})^{1} &\cdots & (w_{n}^{-(n-1)})^{n-1} \end{bmatrix}\begin{bmatrix} A(w_{n}^{0})\\ A(w_{n}^{1})\\ \vdots \\ A(w_{n}^{n-1}) \end{bmatrix}=\begin{bmatrix} a_{0}\\ a_{1}\\ \vdots \\ a_{n-1} \end{bmatrix}$

这就是相当于把 $w_{n}^{i}$ 变成 $w_{n}^{-i}$ 再做一次DFT,然后每个结果再除于n

给出FFT的模板

const int MAXN = 2e5+5;
struct complex {
    double a, b;
    complex(double aa = 0.0, double bb = 0.0) { a = aa; b = bb; }
    complex operator +(const complex &e) { return complex(a + e.a, b + e.b); }
    complex operator -(const complex &e) { return complex(a - e.a, b - e.b); }
    complex operator *(const complex &e) { return complex(a * e.a - b * e.b, a * e.b + b * e.a); }
}x1[MAXN], x2[MAXN], x[MAXN];
void change(complex *y, int len){
    int i, j, k;
    for(i=1,j=len/2; i<len-1; i++){
        if(i < j) swap(y[i], y[j]);
        k = len/2;
        while(j >= k){ j-=k; k>>=1; }
        if(j < k) j += k;
    }
}
void fft(complex *y, int len, int on){
    change(y, len);
    for(int h=2; h<=len; h<<=1){
        complex wn(cos(-on*2*PI/h), sin(-on*2*PI/h));
        for(int j=0; j<len; j+=h){
            complex w(1, 0);
            for(int k=j; k<j+h/2; k++){
                complex u = y[k], t = w*y[k+h/2];
                y[k] = u + t, y[k+h/2] = u- t, w = w*wn;
            }
        }
    }
    if(on == -1) for(int i=0; i<len; i++)  y[i].a /= len;
}

FFT可以来求解大数相乘的问题,一个整数可以拆成n位多项式相加的形式就可以用FFT来求出相乘以后的系数了

如HDU 1402
求两个整数的乘积，两个数的长度有50000位

思路

一个n位数可以拆成 $a_{0}+a_{1}x^1+a_{2}x^{2}+\cdots+a_{n}x^n$ 的形式
x就可以认为是10,那么两个数相乘就符合卷积的运算所以可以用FFT进行运算

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <math.h>
#define PI acos(-1)
using namespace std;
typedef long long ll;
const int MAXN = 2e5+5;
struct complex {
    double a, b;
    complex(double aa = 0.0, double bb = 0.0) { a = aa; b = bb; }
    complex operator +(const complex &e) { return complex(a + e.a, b + e.b); }
    complex operator -(const complex &e) { return complex(a - e.a, b - e.b); }
    complex operator *(const complex &e) { return complex(a * e.a - b * e.b, a * e.b + b * e.a); }
}x1[MAXN], x2[MAXN], x[MAXN];
void change(complex *y, int len){
    int i, j, k;
    for(i=1,j=len/2; i<len-1; i++){
        if(i < j) swap(y[i], y[j]);
        k = len/2;
        while(j >= k){ j-=k; k>>=1; }
        if(j < k) j += k;
    }
}
void fft(complex *y, int len, int on){
    change(y, len);
    for(int h=2; h<=len; h<<=1){
        complex wn(cos(-on*2*PI/h), sin(-on*2*PI/h));
        for(int j=0; j<len; j+=h){
            complex w(1, 0);
            for(int k=j; k<j+h/2; k++){
                complex u = y[k], t = w*y[k+h/2];
                y[k] = u + t, y[k+h/2] = u- t, w = w*wn;
            }
        }
    }
    if(on == -1) for(int i=0; i<len; i++)  y[i].a /= len;
}
char str1[MAXN],str2[MAXN];
long long ans[MAXN];
int main()
{
    while(scanf("%s%s",str1,str2)!=EOF)
    {
        memset(x1,0,sizeof(x1));
        memset(x2,0,sizeof(x2));
        memset(x,0,sizeof(x));
        int len1=strlen(str1);
        int len2=strlen(str2);
        for(int i=0;i<len1;i++)
            {
                x1[len1-1-i].a=str1[i]-'0';
                x1[len1-1-i].b=0;
            }
        for(int i=0;i<len2;i++)
            {
                x2[len2-1-i].a=str2[i]-'0';
                x2[len2-1-i].b=0;
            }
        int len=1;
        while((len<<=1)<(len1+len2));
        fft(x1,len,1);
        fft(x2,len,1);
        for(int i=0;i<len;i++)
            x[i]=x1[i]*x2[i];
        fft(x,len,-1);
        memset(ans,0,sizeof(ans));
        for(int i=0;i<len;i++)
        ans[i]=(ll)(x[i].a+0.5);
        for(int i=0;i<len;i++)
        {
            ans[i+1]+=ans[i]/10;
            ans[i]%=10;
        }
        int temp=0;
        for(int i=len-1;i>=0;i--)
            if(ans[i]!=0)
            {
                temp=i;
                break;
            }
        for(int i=temp;i>=0;i--)
            putchar(ans[i]+'0');
        puts("");
    }
    return 0;
}

兔蠢蠢

关注

1
点赞
踩
8

收藏

觉得还不错? 一键收藏
0
评论
FFT求解卷积

卷积已知两个多项式 A(x)=a0+a1x1+a2x2+a3x3+⋯+anxn" role="presentation" style="position: relative;">A(x)=a0+a1x1+a2x2+a3x3+⋯+anxnA(x)=a0+a1x1+a2x2+a3x3+⋯+anxnA(x)=a_{0}+a_{1}x^1+a_{2}x^2+a_{3}x^3+\cdot
复制链接

扫一扫

专栏目录