当前位置: 首页>后端>正文

python 二维傅里叶逆变换平移 11+w2傅里叶逆变换

  • 前言
  • 三角函数
  • 前前言
  • 单位圆
  • 公式
  • 虚数
  • 定义
  • 公式
  • 单位复根
  • 坐标轴
  • 定义
  • 消去定理
  • 消去定理的推论
  • 折半定理
  • 求和定理
  • 多项式的表达
  • 点值表达
  • 系数表达
  • 系数表达与点值表达的关系
  • 快速傅里叶变换&快速傅里叶逆变换
  • 快速傅里叶变换&快速傅里叶逆变换的优化
  • 最终代码

前言

这篇文章咕的很久,三角函数似乎没啥用。

三角函数

前前言

三角函数似乎没啥用。三角函数似乎没啥用。三角函数似乎没啥用。

单位圆

一个以原点为中心且半径为 \(1\)

公式

\[\cos(\alpha-\beta)=\cos\alpha\cos\beta+\sin\alpha\sin\beta \]

\[\cos(\alpha+\beta)=\cos\alpha\cos\beta-\sin\alpha\sin\beta \]

\[\sin(\alpha+\beta)=\sin\alpha\cos\beta+\cos\alpha\sin\beta \]

\[\sin(\alpha-\beta)=\sin\alpha\cos\beta-\cos\alpha\sin\beta \]

作以 \(x\) 正半轴为始边,逆时针旋转 \(\alpha\),\(\beta\),\(\alpha-\beta\) 的终边与单位圆分别交于 \(A(\cos(\alpha),\sin(\alpha))\),\(B(\cos(\beta),\sin(\beta))\),\(C(\cos(\alpha-\beta),\sin(\alpha-\beta))\), \(x\) 正半轴与单位圆交于 \(D(1,0)\),有 \(\lvert AB \rvert = \lvert CD \rvert\)。

\[(\cos\alpha-\cos\beta)^2+(\sin\alpha-\sin\beta)^2=(1-\cos(\alpha-\beta))^2+\sin^2(\alpha-\beta) \]

\[\begin{aligned} (\cos\alpha-\cos\beta)^2+(\sin\alpha-\sin\beta)^2&=\cos^2\alpha -2\cos\alpha\cos\beta+\cos^2\beta +\sin^2\alpha-2\sin\alpha\sin\beta+\sin^2\beta\\ &=2-2\cos\alpha\cos\beta-2\sin\alpha\sin\beta\\ (1-\cos(\alpha-\beta))^2+\sin^2(\alpha-\beta)&=1-2\cos(\alpha-\beta)+\cos^2(\alpha-\beta)+\sin^2(\alpha-\beta)\\ &=2-2\cos(\alpha-\beta) \end{aligned}\]

\[2-2\cos\alpha\cos\beta-2\sin\alpha\sin\beta=2-2\cos(\alpha-\beta) \]

\[\cos\alpha\cos\beta+\sin\alpha\sin\beta=\cos(\alpha-\beta) \]

\[\cos(-\beta)=\cos\beta \]

\[\sin(-\beta)=-\sin\beta \]

\[\cos(\alpha+\beta)=\cos\alpha\cos\beta-\sin\alpha\sin\beta \]

\[\begin{aligned} \sin(\alpha+\beta)&=\cos(\dfrac{\pi}{2}-\alpha-\beta)=\cos((\dfrac{\pi}{2}-\alpha)-\beta)\\ &=\cos(\dfrac{\pi}{2}-\alpha)\cos\beta+\sin(\dfrac{\pi}{2}-\alpha)\sin\beta\\ &=\sin\alpha\cos\beta+\cos\alpha\sin\beta \end{aligned}\]

\[\begin{aligned} \sin(\alpha-\beta)&=\sin\alpha\cos(-\beta)-\cos\alpha\sin(-\beta)\\ &=\sin\alpha\cos\beta-\cos\alpha\sin\beta \end{aligned}\]

虚数

定义

\[i=\sqrt{-1} \]

\[a+bi \]

公式

\[(a+bi)+(c+di)=(a+c)+(b+d)i \]

\[(a+bi)-(c+di)=(a-c)+(b-d)i \]

\[(a+bi)\times(c+di)=ac+adi+bci+bd\times(\sqrt{-1})^2=(ac-bd)+(ad+bc)i \]

单位复根

坐标轴

\(x\) 代表实数 \(y\)

定义

\(n\) 次单位复根是满足 \(\omega_n^n=1\) 的复数 \(\omega_n\),一共有 \(n\)

如何理解?画一个单位圆,将其 \(n\) 等分,圆上 \(n\) 个点就是 \(n\)

因此有:

\[\omega_n=\cos(\dfrac{2\pi}{n})+\sin(\dfrac{2\pi}{n})i \]

\[\omega_n^k=\cos(\dfrac{2\pi}{n}k)+\sin(\dfrac{2\pi}{n}k)i \]

消去定理

\[\omega_{nm}^{km}=\omega_n^k \]

\[\omega_{nm}^{km}=\cos(\dfrac{2\pi}{nm}km)+\sin(\dfrac{2\pi}{nm}km)i=\cos(\dfrac{2\pi}{n}k)+\sin(\dfrac{2\pi}{n}k)i=\omega_n^k \]

消去定理的推论

\[\omega_n^{n/2}=\omega_2=-1 \]

折半定理

若 'n&1=0' 且 \(n>0\),有:

\[(\omega_n^{k+n/2})^2=(\omega_n^{k})^2 \]

\[(\omega_n^{k+n/2})^2=(\omega_n^{k}\omega_n^{n/2})^2=(\omega_n{^k}\times(-1))^2=(\omega_n^{k})^2 \]

求和定理

对于任意整数 \(1\le n\) 和不能被 \(n\) 整除的非负整数 \(k\),有

\[\sum_{j=0}^{n-1}(\omega_{n}^k)^j=0 \]

多项式的表达

点值表达

一个次数界为 \(n\) 的多项式 \(A(x)\) 的点值表达是一个由 \(n\)

\[\left\{(x_0,y_0),(x_1,y_1),\cdots,(x_{n-2}, y_{n-2}),(x_{n-1},y_{n-1})\right\} \]

其中 \(x_i\)

一个傻逼的性质, \(A(x) = B(x) + C(x)\)

\[A(x) = \left\{(x_0,By_0+Cy_0),(x_1,By_1+Cy_1),\cdots,(x_{n-2},B y_{n-2}+C y_{n-2}),(x_{n-1},By_{n-1}+Cy_{n-1})\right\} \]

\(A(x) = B(x) \times C(x)\)

\[A(x) = \left\{(x_0,By_0\times Cy_0),(x_1,By_1\times Cy_1),\cdots,(x_{n-2},By_{n-2}\times C y_{n-2}),(x_{n-1},By_{n-1}\times Cy_{n-1})\right\} \]

系数表达

一个次数界为 \(n\) 的多项式 \(A(x)\) 的系数表达是一个由 \(n\)

\[a_0x^0 + a_1x^1 + a_2x^2 + \cdots + a_{n-2}x^{n-2}+a_{n-1}x^{n-1} \]

系数表达与点值表达的关系

对于任意一个由 \(n\) 个点对组成的点值表达都有唯一一个次数界为 \(n\)

有 \(y_i = A(x_i)\)。

把其写成矩阵形式:
\( \begin{bmatrix} 1 \quad x_0 \quad x_0^2 \quad &\cdots \quad x_0^{n-1}\\1 \quad x_1 \quad x_1^2 \quad &\cdots \quad x_1^{n-1}\\\vdots\quad\vdots\quad\vdots\quad&\ddots\quad\vdots\\1 \quad x_{n-1} \quad x_{n-1}^2 \quad &\cdots \quad x_{n-1}^{n-1}\\\end{bmatrix}\) \(\begin{bmatrix}a_0\\a_1\\\vdots\\a_{n-1}\end{bmatrix}\) \(=\) \(\begin{bmatrix}y_0\\y_1\\\vdots\\y_{n-1}\end{bmatrix}\)
范德蒙德矩阵在 \(x_i\) 皆不同的情况下有逆矩阵,因此有点值表达,能够确定系数 \(a_i\)。

快速傅里叶变换&快速傅里叶逆变换

规定 \(n\) 是 \(2\) 的 \(q\)

已知

\[A(x)=a_0x^0 + a_1x^1 + a_2x^2 + \cdots + a_{n-2}x^{n-2}+a_{n-1}x^{n-1} \]

\[\left\{(\omega_n,y_0),(\omega_n^1,y_1),\cdots,(\omega_n^{n-2}, y_{n-2}),(\omega_n^{n-1},y_{n-1})\right\} \]

分治!!!

\[A_0(x)=a_0x^0 + a_2x^1 + \cdots + a_{n-2}x^{n/2-1} \]

\[A_1(x)=a_1x^0 + a_3x^1 + \cdots +a_{n-1}x^{n/2-1} \]

\[A_2(x)=x \]

\[A(x)=A_0(x^2)+A_2(x)A_1(x^2) \]

\[A_0=\left\{(\omega_{n/2},y_0),(\omega_{n/2}^1,y_1),\cdots,(\omega_{n/2}^{n/2-2}, y_{n/2-2}),(\omega_{n/2}^{n/2-1},y_{n/2-1})\right\} \]

\[A_1=\left\{(\omega_{n/2},y_0),(\omega_{n/2}^1,y_1),\cdots,(\omega_{n/2}^{n/2-2}, y_{n/2-2}),(\omega_{n/2}^{n/2-1},y_{n/2-1})\right\} \]

\[A_2=\left\{(\omega_{n},y_0),(\omega_{n}^1,y_1),\cdots,(\omega_{n}^{n-2}, y_{n-2}),(\omega_{n}^{n-1},y_{n-1})\right\} \]

\[A(\omega_{n}^{i})=A_0(\omega_{n}^{2i})+\omega_{n}^{i}A_1(\omega_{n}^{2i}) \]

\[A(\omega_{n}^{i})=A_0(\omega_{n/2}^{i})+\omega_{n}^{i}A_1(\omega_{n/2}^{i}) \]

注意,\(\omega_{n/2}^{i}\) 在 \(A_0\) 和 \(A_1\)

巨丑的伪代码:

FFT(a):
    n=a.length
    if n==1
        return a
    w=e^{2 * pi * i / n}
    W=1
    A0={a0,a2,...,a(n-2)}
    A1={a1,a3,...,a(n-1)}
    y0=FFT(A0)
    y1=FFT(A1)
    for k in range(n/2):   // 0 至 n/2-1
        y[k]=y0[k]+W*y1[k]
        y[k]=y0[k]-W*y1[k]
        W*=w
    return y

这样写为什么是对的???

\[\begin{aligned} y0_k&=A0(\omega_{n/2}^{k})\\ y1_k&=A1(\omega_{n/2}^{k})\\ y_k&=A(\omega_{n}^{k})\\ A(\omega_{n}^{k})&=y0_k+\omega_{n}^{k}y1_k\\ &=A0(\omega_{n/2}^{k})+\omega_{n}^{k}A1(\omega_{n/2}^{k})\\ A(\omega_{n}^{k+n/2})&=A0((\omega_{n}^{k+n/2})^2)+\omega_{n}^{k+n/2}A1((\omega_{n}^{k+n/2})^2)\\ &=A0(\omega_{n}^{2k})+\omega_{n}^{n/2}\omega_{n}^{k}A1(\omega_{n}^{2k})\\ &=A0(\omega_{n/2}^{k})+(-1)\omega_{n}^{k}A1(\omega_{n/2}^{k})\\ &=A0(\omega_{n/2}^{k})-\omega_{n}^{k}A1(\omega_{n/2}^{k})\\ \end{aligned} \]

已知

\[\left\{(\omega_n,y_0),(\omega_n^1,y_1),\cdots,(\omega_n^{n-2}, y_{n-2}),(\omega_n^{n-1},y_{n-1})\right\} \]

\[A(x)=a_0x^0 + a_1x^1 + a_2x^2 + \cdots + a_{n-2}x^{n-2}+a_{n-1}x^{n-1} \]

\[ \begin{bmatrix} 1 &1 &1 &1 &1 &\cdots &1\\ 1 &\omega_{n}^{1} &\omega_{n}^{2} &\omega_{n}^{3} &\omega_{n}^{4} &\cdots &\omega_{n}^{n-1}\\ 1 &\omega_{n}^{2} &\omega_{n}^{4} &\omega_{n}^{6} &\omega_{n}^{8} &\cdots &\omega_{n}^{2(n-1)}\\ 1 &\omega_{n}^{3} &\omega_{n}^{6} &\omega_{n}^{9} &\omega_{n}^{12} &\cdots &\omega_{n}^{3(n-1)}\\ 1 &\vdots &\vdots &\vdots &\vdots &\ddots &\vdots\\ 1 &\omega_{n}^{n-1}&\omega_{n}^{2(n-1)}&\omega_{n}^{3(n-1)}&\omega_{n}^{4(n-1)}&\cdots &\omega_{n}^{(n-1)(n-1)} \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ a_2\\ a_3\\ \vdots\\ a_{n-1} \end{bmatrix} = \begin{bmatrix} y_0\\ y_1\\ y_2\\ y_3\\ \vdots\\ y_{n-1} \end{bmatrix} \]

\[V_n= \begin{bmatrix} 1 &1 &1 &1 &1 &\cdots &1\\ 1 &\omega_{n}^{1} &\omega_{n}^{2} &\omega_{n}^{3} &\omega_{n}^{4} &\cdots &\omega_{n}^{n-1}\\ 1 &\omega_{n}^{2} &\omega_{n}^{4} &\omega_{n}^{6} &\omega_{n}^{8} &\cdots &\omega_{n}^{2(n-1)}\\ 1 &\omega_{n}^{3} &\omega_{n}^{6} &\omega_{n}^{9} &\omega_{n}^{12} &\cdots &\omega_{n}^{3(n-1)}\\ 1 &\vdots &\vdots &\vdots &\vdots &\ddots &\vdots\\ 1 &\omega_{n}^{n-1}&\omega_{n}^{2(n-1)}&\omega_{n}^{3(n-1)}&\omega_{n}^{4(n-1)}&\cdots &\omega_{n}^{(n-1)(n-1)} \end{bmatrix}\]

\[V_n^{-1}= \dfrac{1}{n} \begin{bmatrix} 1 &1 &1 &1 &1 &\cdots &1\\ 1 &\omega_{n}^{-1} &\omega_{n}^{-2} &\omega_{n}^{-3} &\omega_{n}^{-4} &\cdots &\omega_{n}^{-(n-1)}\\ 1 &\omega_{n}^{-2} &\omega_{n}^{-4} &\omega_{n}^{-6} &\omega_{n}^{-8} &\cdots &\omega_{n}^{-2(n-1)}\\ 1 &\omega_{n}^{-3} &\omega_{n}^{-6} &\omega_{n}^{-9} &\omega_{n}^{-12} &\cdots &\omega_{n}^{-3(n-1)}\\ 1 &\vdots &\vdots &\vdots &\vdots &\ddots &\vdots\\ 1 &\omega_{n}^{-(n-1)}&\omega_{n}^{-2(n-1)}&\omega_{n}^{-3(n-1)}&\omega_{n}^{-4(n-1)}&\cdots &\omega_{n}^{-(n-1)(n-1)} \end{bmatrix}\]

\[V_n^{-1}V_n=\begin{bmatrix} 1 &0 &0 &\cdots &0\\ 0 &1 &0 &\cdots &0\\ 0 &0 &1 &\cdots &0\\ 0 &\vdots &\vdots &\ddots &\vdots\\ 0 &0 &0 &\cdots &1 \end{bmatrix}\]

\[[V_n^{-1}V_n]_{i,j}=\sum_{k=0}^{n-1} (\omega_n^{-ki}/n)(\omega_n^{kj}) =\sum_{k=0}^{n-1}\omega_n^{kj-ki}/n \]

应用求和引理。

在看原来的柿子:

\[\left\{(\omega_n,y_0),(\omega_n^1,y_1),\cdots,(\omega_n^{n-2}, y_{n-2}),(\omega_n^{n-1},y_{n-1})\right\} \]

\[A(x)=a_0x^0 + a_1x^1 + a_2x^2 + \cdots + a_{n-2}x^{n-2}+a_{n-1}x^{n-1} \]

\[a_j=\dfrac{1}{n}\sum_{k=0}^{n-1}y_k\omega_{n}^{-kj} \]

稍微变换一下 \(FFT\)

快速傅里叶变换&快速傅里叶逆变换的优化

改成迭代计算即可,因为在同一层的计算会用到相同的单位复根,可以减少常数。

简单的说,将序列位逆序置换,得到最低层的计算顺序,然后不断向上计算,调整位置(实际上在合并时就自动调整位置了)。

最终代码

#include <bits/stdc++.h>
using namespace std;
const int N = 2e7+11;
const double pi = acos(-1);
struct Comp{
    double a, b;
    Comp operator + (const Comp &oth)const{return (Comp){a+oth.a, b+oth.b};};
    Comp operator - (const Comp &oth)const{return (Comp){a-oth.a, b-oth.b};};
    Comp operator * (const Comp &oth)const{return (Comp){a*oth.a-b*oth.b, a*oth.b+b*oth.a};}
}f[N], g[N], w_i[N];
int n, m, rev[N];
void FFT(int len, Comp *a, bool type){
    for(int i = 0; i < len; i ++){
        rev[i] = (rev[i>>1]>>1) + (i%2?(len>>1):0);
        if(rev[i] > i) swap(a[rev[i]], a[i]);
    }
    for(int d = 1; d < len; d<<=1){
        Comp w_n = (Comp){cos(pi/d), sin(pi/d)};
        if(type) w_n.b = -w_n.b;
        w_i[0] = (Comp){1, 0};
        for(int i = 1; i < d; i ++) w_i[i] = w_i[i-1]*w_n;
        for(int fir = 0; fir < len; fir += d<<1){
            int sec = fir+d;
            for(int i = 0; i < d; i ++){
                Comp A0 = a[fir+i], A1 = w_i[i] * a[sec+i];
                a[fir+i] = A0 + A1, a[sec+i] = A0 - A1;
            }
        }
    }
    if(type) for(int i = 0; i < len; i ++)
        a[i].a /= len, a[i].b /= len;
    return;
}
int main(){
    scanf("%d %d", &n, &m);
    for(int i = 0; i <= n; i ++) scanf("%lf", &f[i].a);
    for(int i = 0; i <= m; i ++) scanf("%lf", &g[i].a);
    int len = 1; while(len <= n+m) len <<= 1; cout << len << endl;
    FFT(len, f, 0), FFT(len, g, 0);
    for(int i = 0; i < len; i ++) f[i] = f[i] * g[i], printf("%.2lf ", f[i].a);
    puts("");
    FFT(len, f, 1);
    for(int i = 0; i <= n+m; i ++)
        printf("%d ", (int)(f[i].a+0.5));
    return 0;
}




https://www.xamrdz.com/backend/3uz1960889.html

相关文章: