多项式

一个次数界为的多项式
次数界为说明这个多项式最高次项的指数小于

系数表示法

我们可以只保留多项式每项的系数来描述这个多项式。即表示指数为的项的系数

运算

多项式之间的加运算,只要将对应项的系数相加即可。
即若那么有
多项式之间的乘法。
那么
其中叫做卷积
可以看出暴力计算卷积的复杂度为

点值表示法

我们可以选个值代入多项式从而得出个二元组。可以证明,通过这个点对可以确定这个多项式。
还原这个多项式用高斯消元即可。
点值表示法的优点是进行多项式的加法和乘法只要将相加或者相乘即可。

求值与插值

考虑如何可以比较快速的计算卷积。因为点值表示计算乘法比较方便。所以我们可以把系数表示法转换为点值表示法。然后进行乘法。然后在转换回来。就可以计算卷积了。
把系数表示法转化为点值表示法称为求值,把点值表示法转换为系数表示法称为插值。

拉格朗日插值公式


利用拉格朗日插值公式可以在给定个点值表示法的情况下时间内计算出原二项式在某个位置的值

离散傅里叶变换(DFT)

离散傅里叶变换实现了两种表示法之间的转换,复杂度是

复数

复数是形如的数,其中为实数。称为实部,称为虚部。为虚数单位,定义
复数是数的范围在实数上的扩展
将复数表示在数轴上,是一个的向量.

复数的运算

复数相加将分别相加即可。
复数相乘,其实将i看作未知量相乘就行了。即

单位圆

在数轴上以原点为圆心半径为的圆。

单位根

次单位根满足
表现在单位圆上就是单位圆的等分点
我们用表示主次单位根也就是从开始数第一个n次单位根
那么可以得到其他的次单位根分别为
性质: 当n为偶数时,

欧拉公式



用来求主次单位根

转换

有了上面这些铺垫,就可以进行转换了。
次单位根带入原来的多项式,就可以得到个点值表示法。
那么怎么将点值表示法转换为系数表示法呢?
这就是用次单位根的原因
我们将得到作为另一个多项式的系数
然后取单位根倒数带入就能得到原来的多项式
这样就成功的实现了点值表示法与系数表示法之间的转换

快速傅里叶变换(FFT)

计算卷积的复杂度是的。利用快速傅里叶变换可以优化到

分治

考虑一个次数界为的多项式求他的
这里默认为偶数,如果不是偶数,那么只要加上一个系数为的项即可将他填充为偶数
然后我们根据下标的奇偶性将他分为两部分。
$A(x)=A_{(even)}(x^2) + xA_{(odd)}(x^2)A_{(even)}(x)A_{(odd)}(x)n\omega_{n}^iA(\omega_n^{i + \frac{n}{2}})n$个单位根代入所得到的值

优化

递归实现非常慢。所以要用非递归的方法来实现
我们将递归实现的过程写下来
0 1 2 3 4 5
0 2 4|1 3 5
0 4|2|1 5|3
0|4|2|1|5|3
然后看出每个数字递归到最后的为是为当前位置二进制表示法翻转之后的数字。
然后就可以先把每个数字都放到最后的位置上面,然后再合并就行了。

代码

递归

/*
* @Author: wxyww
* @Date:   2019-02-01 21:09:51
* @Last Modified time: 2019-02-01 21:47:34
*/
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<cmath>
#include<ctime>
#include<cmath>
#include<bitset>
using namespace std;
typedef long long ll;
const int N = 3000100;
const double pi = acos(-1.0);
ll read() {
    ll x=0,f=1;char c=getchar();
    while(c<'0'||c>'9') {
        if(c=='-') f=-1;
        c=getchar();
    }
    while(c>='0'&&c<='9') {
        x=x*10+c-'0';
        c=getchar();
    }
    return x*f;
}
struct complex {
    double x,y;
    complex() {x = y = 0;}
    complex(double xx,double yy) {
        x = xx;y = yy;
    }
}A[N],B[N];
complex operator * (complex a,complex b) {
    return complex(a.x * b.x - a.y * b.y,a.x * b.y + a.y * b.x);//复数运算
}
complex operator + (complex a,complex b) {
    return complex(a.x + b.x,a.y + b.y);
}
complex operator - (complex a,complex b) {
    return complex(a.x - b.x,a.y - b.y);
}

void FFT(complex *a,int n,int ty) {
    if(n == 1) return;
    complex a1[n >> 1],a2[n >> 1];
    for(int i = 0;i <= n;i += 2) {
        a1[i >> 1] = a[i];a2[i >> 1] = a[i + 1];//按奇偶分开
    }
    FFT(a1,n >> 1,ty);FFT(a2,n >> 1,ty);//递归
    complex w1 = complex(cos(2.0 * pi / n),ty * sin(2.0 * pi / n));//主n次单位根
    complex w = complex(1.0,0.0);//当前的n次单位根
    int k = n >> 1;
    for(int i = 0;i < k;++i) {//根据式子计算
        complex t = w * a2[i];
        a[i + k] = a1[i] - t;
        a[i] = a1[i] + t;
        w = w *  w1;
    }

}
int main() {
    int n = read(),m = read();
    for(int i = 0;i <= n;++i) A[i].x = read();
    for(int i = 0;i <= m;++i) B[i].x = read();
    int tot = 1;
    while(tot <= n + m) tot <<= 1;
    FFT(A,tot,1);
    FFT(B,tot,1);
    for(int i = 0;i <= tot;++i) A[i] = A[i] * B[i];
    FFT(A,tot,-1);
    for(int i = 0;i <= n + m;++i) 
        printf("%d ",(int)(A[i].x / tot + 0.5));
    return 0;
}

非递归

/*
* @Author: wxyww
* @Date:   2019-02-01 21:50:51
* @Last Modified time: 2019-02-01 22:03:10
*/
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<cmath>
#include<ctime>
#include<cmath>
#include<bitset>
using namespace std;
typedef long long ll;
const int N = 3000100;
const double pi = acos(-1.0);
ll read() {
    ll x=0,f=1;char c=getchar();
    while(c<'0'||c>'9') {
        if(c=='-') f=-1;
        c=getchar();
    }
    while(c>='0'&&c<='9') {
        x=x*10+c-'0';
        c=getchar();
    }
    return x*f;
}
struct complex {
    double x,y;
    complex() {x = y = 0;}
    complex(double xx,double yy) {
        x = xx;y = yy;
    }
}A[N],B[N];
complex operator * (complex a,complex b) {
    return complex(a.x * b.x - a.y * b.y,a.x * b.y + a.y * b.x);//复数运算
}
complex operator + (complex a,complex b) {
    return complex(a.x + b.x,a.y + b.y);
}
complex operator - (complex a,complex b) {
    return complex(a.x - b.x,a.y - b.y);
}

void FFT(complex *a,int n,int ty) {
    for(int i = 0,j = 0;i < n;++i) {//找到最终对应的位置
        if(i < j) swap(a[i],a[j]);
        for(int k = n >> 1;(j ^= k) < k;k >>= 1);
    }
    for(int m = 2;m <= n;m <<= 1) {
        complex w1 = complex(cos(2*pi/m),ty * sin(2 * pi / m));
        for(int i = 0;i < n;i += m) {
            complex w = complex(1,0);
            for(int k = 0;k < (m >> 1);++k) {
                complex t = w * a[i + k + (m >> 1)];
                complex u = a[i + k];
                a[i + k] = u + t;
                a[i + k + (m >> 1)] = u - t;
                w = w * w1;
            }
        }
    }
}
int main() {
    int n = read(),m = read();
    for(int i = 0;i <= n;++i) A[i].x = read();
    for(int i = 0;i <= m;++i) B[i].x = read();
    int tot = 1;
    while(tot <= n + m) tot <<= 1;
    FFT(A,tot,1);
    FFT(B,tot,1);
    for(int i = 0;i <= tot;++i) A[i] = A[i] * B[i];
    FFT(A,tot,-1);
    for(int i = 0;i <= n + m;++i) 
        printf("%d ",(int)(A[i].x / tot + 0.5));
    return 0;
}

时间差异

最后来一份递归代码与非递归代码运行时间的比较