FFT多项式乘法学习笔记

??其实我不知道我是否真的理解了FFT,但是我会用FFT优化多项式乘法了QAQ。。

(以下大多摘自算导

前置知识

1. 多项式

??在一个代数域F上,关于变量x的多项式定义为形式和形式表示的函数

A(x)=∑j=0n?1ajxj,其中a0…an?1为多项式各项的系数

2. 多项式的次数界

??若多项式有非零系数的最高次项为xk,则称k为该多项式的次数,任何严格大于k的整数都是这个多项式的次数界。

3. 多项式的表示

(1)系数表示法

??对于一个次数界为n的多项式A(x)来说,其系数表示法可以看做是一个列向量a=(a0,a1,…,an?1)。系数表示法对于某些多项式的计算很方便,如对多项式A(x)在给定点x0的求值运算就是计算A(x0)的值。如果使用霍纳法则,则求值运算的运行时间为O(n),

A(x0)=a0+x0(a1+x0(a2+...+x0(an?2+x0(an?1))...))

??另外,加法运算的时间复杂度是O(n),暴力进行乘法运算的时间复杂度是O(n2)。

(2)点值表示法

?? 有n个点(x0,A(x0),),(x1,A(x1)),..,(xn?1,A(n?1)),当所有xk各不相同时,这n个点可以唯一表示一个次数界为n的多项式,但是一个次数界为n的多项式可以有多个点值表示。通俗一点说,已知一个多项式函数的n个函数值可以唯一确定这个函数,而知道这个函数可以知道不止n个函数值。

??已知点值表示求系数表示称为插值,用拉格朗日插值法可以做到O(n2)。

??若次数界为n的多项式A(x),B(x)的n个点的xk是对应相同的,点值表示法的加法操作时间复杂度是O(n),只要把对应A/B(xk)相加即可,若A(x),B(x)都已知2n个点,那把A/B(xk)相乘同样可以在O(n)时间内完成乘法运算。。

进入正题。。

1. 单位复根

??n次单位复根是满足wn=1的复数w,有n个,他们均匀的分布在以复

平面的原点为圆心的单位圆上,为e2πki/n,k=0,1,…,n?1,复数幂定义为eui=cos(u)+sin(u)i 。wn=e2πi/n称为主n次单位根,所有的其他n次单位根都是wn的幂。

??因为有wnn=w0n=1 ,所以有wjnwkn=w(j+k)modnn 。(性质1

重要性质

相消引理

?? 对任何整数n≥0,k≥0,d>0,有 wdkdn=wkn,利用定义不难证明

?? 推论:对任意偶数n>0,有wn/2n=w2=1 。

折半引理

?? 如果n>0为偶数,n个n次单位复根的平方等于n/2个n/2次单位复根。利用相消引理和性质1可证。

求和引理

??对于任意整数n≥1和不能被n整除的非零整数k,有 ∑n?1j=0wkjn=0。

证明:

?? 原式=(wkn)n?1wkn?1=(wnn)k?1wkn?1=1?1wkn?1=0 ,只要k不整除n就可以保证分母不为0。

2. DFT

??A(x)是一个次数界为n的多项式,不失一般性地假定n是2的幂,因为

次数界总是可以增大的。有列向量y=(y0,y1,…,yn?1),其中yk=A(wkn),则称y是系数向量a的离散傅里叶变换,也写作y=DFTn(a)。

3. 快速傅里叶变换(FFT)

??FFT是一种可以在O(nlogn)时间内计算出DFTn(a)的算法,主要思想

是分治。

??定义A0(x)=a0+a2x+a4x2+…+an?2xn/2?1,包含了A(x)

所有偶数下标的系数,A1(x)=a1+a3x+a5x2+…+an?1xn/2?1,包含了A(x)所有奇数下标的系数。易得A(x)=A0(x2)+xA1(x2),所以我们可以先求A0和A1的DFT,然后再组合起来。组合的时候有

A(wk+n/2n)=A0(w2k+nn)+wk+n/2nA1(w2k+nn)=A0(wkn)+wknw2A1(wkn)A0(wkn)?wknA1(wkn)

所以用n2个yk就可以推出全部。根据折半引理,问题的规模缩小一半,每次组合的时间复杂度是O(n),所以总时间复杂度是O(nlogn)。

?? 求出DFTn(a)后,要对单位复根进行插值,将点值表示转化为系数表示,有y=Vna,其中

Vn=????????????1111...11wnw2nw3n...wn?1n1w2nw4nw6n...w2n1w3nw6nw9n...w3n..................1wn?1nw2(n?1)nw3(n?1)n...w(n?1)(n?1)n????????????

是范德蒙特矩阵。

a=Vn?1?y,考虑Vn?1在(j,k)处的数Vj,k,根据V?1nVn=In(n阶单位矩阵)有

∑t=0n?1Vj,twktn={1(j=k)0(j!=k)

定理:Vj,k=w?kjn/n

证明:∑n?1t=0Vj,twktn=1n∑n?1t=0wt(k?j)n,当j=k时,上式=1,否则由求和引理可得上式为0,得证。

??那么有 aj=1n∑n?1k=0ykw?kjn,因此在求逆DFT的时候可以类似于求DFT,只要把ya角色互换,然后让 wn=w?1n,再做FFT即可。

??写的时候发现非递归要比递归快很多。。

递归:

#include<cstdio>
#include<iostream>
#include<cmath>
#include<memory.h>
#define N 400010
using namespace std;
const double pi=acos(-1);
struct complex{
    double x,i;
    complex(){}
    complex(double x,double i):x(x),i(i){}
    complex operator+(complex a) {return complex(x+a.x,i+a.i);}
    complex operator-(complex a) {return complex(x-a.x,i-a.i);}
    complex operator*(complex a) {return complex(x*a.x-i*a.i,x*a.i+i*a.x);}
}a[N],b[N];
int n,m,i,nn;
void fft(complex *a,int n,int t)
{
    if (n==1) return;
    complex a0[n>>1],a1[n>>1];
    for (int i=0;i<n;i+=2) a0[i>>1]=a[i],a1[i>>1]=a[i+1];
    fft(a0,n>>1,t);fft(a1,n>>1,t);
    complex wn(cos(2*pi/n),t*sin(2*pi/n)),w(1,0);
    for (int i=0;i<(n>>1);i++,w=w*wn) a[i]=a0[i]+w*a1[i],a[i+(n>>1)]=a0[i]-w*a1[i];
}
int main()
{
    scanf("%d%d",&n,&m);
    memset(a,0,sizeof(a));memset(b,0,sizeof(b));
    for (i=0;i<=n;i++) scanf("%lf",&a[i].x);
    for (i=0;i<=m;i++) scanf("%lf",&b[i].x);
    nn=1;while (nn<=n+m) nn<<=1;
    fft(a,nn,1);fft(b,nn,1);
    for (i=0;i<=nn;i++) a[i]=a[i]*b[i];
    fft(a,nn,-1);
    for (i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/nn+0.5));
}

非递归:

#include<cstdio>
#include<iostream>
#include<cmath>
#include<memory.h>
#define N 400010
using namespace std;
const double pi=acos(-1);
struct complex{
    double x,i;
    complex(){}
    complex(double x,double i):x(x),i(i){}
    complex operator+(complex a) {return complex(x+a.x,i+a.i);}
    complex operator-(complex a) {return complex(x-a.x,i-a.i);}
    complex operator*(complex a) {return complex(x*a.x-i*a.i,x*a.i+i*a.x);}
}a[N],b[N];
int n,m,i,nn,len,rev[N];
void fft(complex *a,int n,int t)
{
    for (int i=0;i<n;i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
    for (int j=1;j<n;j<<=1)
    {
        complex wn(cos(2*pi/(j<<1)),t*sin(2*pi/(j<<1)));
        for (int i=0;i<n;i+=(j<<1))
        {
            complex w(1,0),t0,t1;
            for (int k=0;k<j;k++,w=w*wn) t0=a[i+k],t1=w*a[i+j+k],a[i+k]=t0+t1,a[i+j+k]=t0-t1;
        }
    }
}
int main()
{
    freopen("FFT.in","r",stdin);
    scanf("%d%d",&n,&m);
    memset(a,0,sizeof(a));memset(b,0,sizeof(b));
    for (i=0;i<=n;i++) scanf("%lf",&a[i].x);
    for (i=0;i<=m;i++) scanf("%lf",&b[i].x);
    nn=1;len=0;while (nn<=n+m) nn<<=1,len++;
    rev[0]=0;
    for (i=1;i<nn;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
    fft(a,nn,1);fft(b,nn,1);
    for (i=0;i<=nn;i++) a[i]=a[i]*b[i];
    fft(a,nn,-1);
    for (i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/nn+0.5));
}

??非递归这里有一个翻转的函数,意在让所有数按合并时候的二叉树的叶子节点的顺序排列,不难发现翻转过来就是它的新位置。。

时间: 2024-08-11 05:35:09

FFT多项式乘法学习笔记的相关文章

FFT多项式乘法加速

FFT基本操作...讲解请自己看大学信号转置系列... 15-5-30更新:改成结构体的,跪烂王学长啊啊啊啊机智的static... 1 #include<iostream> 2 #include<cstdio> 3 #include<cmath> 4 #include<algorithm> 5 #include<queue> 6 #include<cstring> 7 #define PAU putchar(' ') 8 #defi

FFT多项式乘法模板

有时间来补算法原理orz #include <iostream> #include <cstdio> #include <cmath> #include <complex> using namespace std; const double pi = acos(-1); const int maxn = 111111; typedef complex<double> Complex; void DFT(Complex *a, int n, int

多项式乘法(FFT)学习笔记

------------------------------------------本文只探讨多项式乘法(FFT)在信息学中的应用如有错误或不明欢迎指出或提问,在此不胜感激 多项式 1.系数表示法     一般应用最广泛的表示方式     用A(x)表示一个x-1次多项式,a[i]为$ x^i$的系数,则A(x)=$ \sum_0^{n-1}$ a[i] * $ x^i$ 仅利用这种方式求多项式乘法复杂度为O($ n^2$),不够优秀2.点值表示法     将n个互不相同的值$ x_0$...$

[学习笔记] 多项式与快速傅里叶变换(FFT)基础

引入 可能有不少OIer都知道FFT这个神奇的算法, 通过一系列玄学的变化就可以在 $O(nlog(n))$ 的总时间复杂度内计算出两个向量的卷积(或者多项式乘法/高精度乘法), 而代码量却非常小. 博主一年半前曾经因COGS的一道叫做"神秘的常数 $\pi$"的题目而去学习过FFT, 但是基本就是照着板子打打完并不知道自己在写些什么鬼畜的东西OwO 不过...博主这几天突然照着算法导论自己看了一遍发现自己似乎突然意识到了什么OwO然后就打了一道板子题还1A了OwO再加上午考试差点AK

【笔记篇】(理论向)快速傅里叶变换(FFT)学习笔记w

现在真是一碰电脑就很颓废啊... 于是早晨把电脑锁上然后在旁边啃了一节课多的算导, 把FFT的基本原理整明白了.. 但是我并不觉得自己能讲明白... Fast Fourier Transformation, 快速傅里叶变换, 是DFT(Discrete Fourier Transform, 离散傅里叶变换)的快速实现版本. 据说在信号处理领域广泛的应用, 而且在OI中也有广泛的应用(比如SDOI2017 R2至少考了两道), 所以有必要学习一波.. 划重点: 其实学习FFT最好的教材是<算法导论

概率图模型(PGM)学习笔记(四)-贝叶斯网络-伯努利贝叶斯-多项式贝叶斯

指针悬空 指针悬空在我们使用指针的时候很容易被忽视,主要的表现是:指针所指向的内存 释放,指针并没有置为NULL,致使一个不可控制的指针. #include<stdio.h> #include<stdlib.h> int *pointer; void func() { int n=8; pointer=&n; printf("pointer point data is %d\n",*pointer); // pointer=NULL; } int mai

Servlet乘法表学习笔记

一.控制台实现乘法表 package com.shanrengo; import java.io.IOException; import java.io.PrintWriter; import javax.servlet.ServletException; import javax.servlet.http.HttpServlet; import javax.servlet.http.HttpServletRequest; import javax.servlet.http.HttpServle

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

P3803 [模板]多项式乘法(FFT) 题目背景 这是一道FFT模板题 题目描述 给定一个n次多项式F(x),和一个m次多项式G(x). 请求出F(x)和G(x)的卷积. 输入输出格式 输入格式: 第一行2个正整数n,m. 接下来一行n+1个数字,从低到高表示F(x)的系数. 接下来一行m+1个数字,从低到高表示G(x))的系数. 输出格式: 一行n+m+1个数字,从低到高表示F(x)∗G(x)的系数. 输入输出样例 输入样例#1: 复制 1 2 1 2 1 2 1 输出样例#1: 复制 1

洛谷.3803.[模板]多项式乘法(FFT)

题目链接:洛谷.LOJ. FFT相关:快速傅里叶变换(FFT)详解.FFT总结.从多项式乘法到快速傅里叶变换. #include <cmath> #include <cctype> #include <cstdio> #include <algorithm> #define gc() getchar() const int N=1e6+5; const double PI=acos(-1); int n,m; struct Complex { double