码迷,mamicode.com
首页 > 其他好文 > 详细

FFT快速傅里叶变换

时间:2020-02-09 20:10:38      阅读:75      评论:0      收藏:0      [点我收藏+]

标签:lin   target   element   show   ide   不同的   omega   ++   minus   

摘自:https://www.cnblogs.com/RabbitHu/p/FFT.html

快速傅里叶变换(FFT)是一种能在O(nlogn)O(nlog?n)的时间内将一个多项式转换成它的点值表示的算法。

点值表示:设A(x)是一个n1次多项式,那么把n个不同的x代入,会得到ny。n对(x,y)唯一确定了该多项。由多项式可以求出其点值表示,而由点值表示也可以求出多项式。

设有两个n1次多项式A(x)B(x)),我们的目标是——把它们乘起来。普通的多项式乘法是O(n^2),但有趣的是,两个用点值表示的多项式相乘,复杂度是O(n)的!具体方法:C(xi)=A(xi)×B(xi)

 

高精度乘法:等价于x=10的A(x)*B(x)

技术图片
 1 #include <cstdio>
 2 #include <cmath>
 3 #include <cstring>
 4 #include <algorithm>
 5 #include <complex>
 6 #define space putchar(‘ ‘)
 7 #define enter putchar(‘\n‘)
 8 using namespace std;
 9 typedef long long ll;
10 template <class T>
11 void read(T &x){
12     char c;
13     bool op = 0;
14     while(c = getchar(), c < 0 || c > 9)
15     if(c == -) op = 1;
16         x = c - 0;
17     while(c = getchar(), c >= 0 && c <= 9)
18         x = x * 10 + c - 0;
19     if(op) x = -x;
20 }
21 template <class T>
22 void write(T x){
23     if(x < 0) putchar(-), x = -x;
24     if(x >= 10) write(x / 10);
25     putchar(0 + x % 10);
26 }
27 const int N = 1000005;
28 const double PI = acos(-1);
29 typedef complex <double> cp;
30 char sa[N], sb[N];
31 int n = 1, lena, lenb, res[N];
32 cp a[N], b[N], omg[N], inv[N];
33 void init(){
34     for(int i = 0; i < n; i++){
35         omg[i] = cp(cos(2 * PI * i / n), sin(2 * PI * i / n));
36         inv[i] = conj(omg[i]);
37     }
38 }
39 void fft(cp *a, cp *omg){
40     int lim = 0;
41     while((1 << lim) < n) lim++;
42     for(int i = 0; i < n; i++){
43         int t = 0;
44         for(int j = 0; j < lim; j++)
45             if((i >> j) & 1) t |= (1 << (lim - j - 1));
46         if(i < t) swap(a[i], a[t]); // i < t 的限制使得每对点只被交换一次(否则交换两次相当于没交换)
47     }
48     for(int l = 2; l <= n; l *= 2){
49         int m = l / 2;
50     for(cp *p = a; p != a + n; p += l)
51         for(int i = 0; i < m; i++){
52             cp t = omg[n / l * i] * p[i + m];
53             p[i + m] = p[i] - t;
54             p[i] += t;
55         }
56     }
57 }
58 int main(){
59     scanf("%s%s", sa, sb);
60     lena = strlen(sa), lenb = strlen(sb);
61     while(n < lena + lenb) n *= 2;
62     for(int i = 0; i < lena; i++)
63         a[i].real(sa[lena - 1 - i] - 0);
64     for(int i = 0; i < lenb; i++)
65         b[i].real(sb[lenb - 1 - i] - 0);
66     init();
67     fft(a, omg);
68     fft(b, omg);
69     for(int i = 0; i < n; i++)
70         a[i] *= b[i];
71     fft(a, inv);
72     for(int i = 0; i < n; i++){
73         res[i] += floor(a[i].real() / n + 0.5);
74         res[i + 1] += res[i] / 10;
75         res[i] %= 10;
76     }
77     for(int i = res[lena + lenb - 1] ? lena + lenb - 1: lena + lenb - 2; i >= 0; i--)
78         putchar(0 + res[i]);
79     enter;
80     return 0;
81 }
82 
83 //可以预处理ωkn和ω−kn,分别存在omg和inv数组中。调用fft时,如果无需取倒数,则传入omg;如果需要取倒数,则传入inv。
View Code

 

FFT快速傅里叶变换

标签:lin   target   element   show   ide   不同的   omega   ++   minus   

原文地址:https://www.cnblogs.com/lbz007oi/p/12288348.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!