快速傅里叶变换

快速傅里叶变换

性质一：可以用n+1个点，表示一个n次多项式

证明用高斯消元，范德蒙行列式满秩唯一解。

点表示法：

那么：如果A(x)是n次的，B(x)是m次的，那么我们能用n+m+1个点表示C(x)。

我们知道C(x)=A(x)B(x)

所以我们可以通过O(n+m)表示出来C(x);

所以我们希望从系数表示法转化成点表示法，在从点表示法转化成系数表示法。

复数

欧拉公式证明：

证明过程参考繁凡博客吧。​​​ using namespace std;const double pi = acos(-1);const int N = 300;struct Complex{ double x, y; Complex(double x = 0, double y = 0): x(x), y(y) {}} A[N], B[N];Complex operator * (Complex J, Complex Q){ //模长相乘，幅度相加 return Complex(J.x * Q.x - J.y * Q.y, J.x * Q.y + J.y * Q.x);}Complex operator - (Complex J, Complex Q){ return Complex(J.x - Q.x, J.y - Q.y);}Complex operator + (Complex J, Complex Q){ return Complex(J.x + Q.x, J.y + Q.y);}void FFT(int limit, Complex *a, int type){ if(limit == 1) return ;//常数项 Complex a1[limit / 2], a2[limit / 2]; //分成两端，左右两端平分 for(int i = 0; i <= limit; i += 2) //偶数项,奇数项 { a1[i / 2] = a[i]; a2[i / 2] = a[i + 1]; } FFT(limit / 2, a1, type); FFT(limit / 2, a2, type); Complex tmp = Complex(cos(2.0 * pi / limit), type * sin(2.0 * pi / limit)), w = Complex(1, 0); ///tmp表示单位根，w表示幂。 for(int i = 0; i < (limit /2); i++, w = w * tmp) { a[i] = a1[i] + w * a2[i]; //左加 a[i + limit / 2] = a1[i] - w * a2[i]; //右减 }}int main(){ int n, m; scanf("%d%d", &n, &m); for(int i = 0; i <= n; i++)scanf("%lf", &A[i].x); for(int i = 0; i <= m; i++)scanf("%lf", &B[i].x); int limit = 1; while(limit <= n + m) limit *= 2; //2的整数次幂 FFT(limit, A, 1); FFT(limit, B, 1); //后面的1表示要进行的变换是什么类型 // 1表示从系数变为点值 //-1表示从点值变为系数 //这就与推导过程中的 w指数部分正负有关。 for(int i = 0; i <= limit; i++) { A[i] = A[i] * B[i]; } FFT(limit, A, -1); for(int i = 0;i <= n + m; i++) { printf("%d ", (int)(A[i].x / limit + 0.5)); ///对其四舍五入+0.5 } return 0;}

迭代版

蝴蝶变换

///迭代版#include using namespace std;const int N = 5e6 + 7;const double PI = acos(-1);int n, m;int limit = 1;int res, ans[N];int l;int r[N];struct Complex{ double x, y; Complex (double x = 0, double y = 0) : x(x), y(y) {}} a[N], b[N];Complex operator * (Complex J, Complex Q){ //模长相乘，幅度相加 return Complex(J.x * Q.x - J.y * Q.y, J.x * Q.y + J.y * Q.x);}Complex operator + (Complex J, Complex Q){ return Complex(J.x + Q.x, J.y + Q.y);}Complex operator - (Complex J, Complex Q){ return Complex(J.x - Q.x, J.y - Q.y);}void FFT(Complex *A, int type){ for(int i = 0; i < limit; i++) { if(i < r[i]) swap(A[i], A[r[i]]); ///保证不在换过去 } ///从底层开始合并： for(int mid = 1; mid < limit; mid = mid * 2) { ///待合并区间长度的一半，最开始是长度为1的合并，类似倍增的思想 Complex wn(cos(PI / mid), type * sin(PI / mid)); //单位根 for(int len = mid *2, pos = 0; pos < limit; pos += len) { ///len是区间的长度，pos是当前的位置 Complex w(1, 0); for(int k = 0; k < mid; k++, w = w * wn) { ///只扫左半部分，蝴蝶变换得到有半部分。 Complex x = A[pos + k]; //左半部分 Complex y = w * A[pos + mid + k]; //右半部分 A[pos + k] = x + y;//左加 A[pos + mid + k] = x - y;///右减 } } } if(type == 1) return ; for(int i = 0; i <= limit; i++) { A[i].x = A[i].x / limit; A[i].y = A[i].y / limit;///这个版本没用，加不加没影响，优化的版本必须要加 ///除以我们推出的N。 }}int main(){ cin >> n >> m; for(int i = 0; i <= n; i++) scanf("%lf", &a[i].x); for(int j = 0; j <= m; j++) scanf("%lf", &b[j].x); while(limit <= n + m) limit <<= 1, l++;///求出位数l，和2的整数次幂 for(int i = 0; i < limit; i++) { r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));///蝴蝶变换。 } FFT(a, 1); FFT(b, 1); for(int i = 0; i <= limit; ++i) { a[i] = a[i] * b[i]; } FFT(a, -1); for(int i = 0; i <= m + n; i++) { printf("%d ", int(a[i].x + 0.5)); } return 0;}

优化：三步变两步：

///优化版#include using namespace std;const int N = 5e6 + 7;const double PI = acos(-1);int n, m;int limit = 1;int res, ans[N];int l;int r[N];struct Complex{ double x, y; Complex (double x = 0, double y = 0) : x(x), y(y) {}} a[N], b[N];Complex operator * (Complex J, Complex Q){ //模长相乘，幅度相加 return Complex(J.x * Q.x - J.y * Q.y, J.x * Q.y + J.y * Q.x);}Complex operator + (Complex J, Complex Q){ return Complex(J.x + Q.x, J.y + Q.y);}Complex operator - (Complex J, Complex Q){ return Complex(J.x - Q.x, J.y - Q.y);}void FFT(Complex *A, int type){ for(int i = 0; i < limit; i++) { if(i < r[i]) swap(A[i], A[r[i]]); ///保证不在换过去 } ///从底层开始合并： for(int mid = 1; mid < limit; mid = mid * 2) { ///待合并区间长度的一半，最开始是长度为1的合并，类似倍增的思想 Complex wn(cos(PI / mid), type * sin(PI / mid)); //单位根 for(int len = mid *2, pos = 0; pos < limit; pos += len) { ///len是区间的长度，pos是当前的位置 Complex w(1, 0); for(int k = 0; k < mid; k++, w = w * wn) { ///只扫左半部分，蝴蝶变换得到有半部分。 Complex x = A[pos + k]; //左半部分 Complex y = w * A[pos + mid + k]; //有半部分 A[pos + k] = x + y; A[pos + mid + k] = x - y; } } } if(type == 1) return ; for(int i = 0; i <= limit; i++) { A[i].x = A[i].x / limit; A[i].y = A[i].y / limit; ///除以我们推出的N。 }}int main(){ cin >> n >> m; for(int i = 0; i <= n; i++) scanf("%lf", &a[i].x); for(int j = 0; j <= m; j++) scanf("%lf", &a[j].y); while(limit <= n + m) limit <<= 1, l++; for(int i = 0; i < limit; i++) { r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); } FFT(a, 1); ///FFT(b, 1); for(int i = 0; i <= limit; ++i) { a[i] = a[i] * a[i]; } FFT(a, -1); for(int i = 0; i <= m + n; i++) { printf("%d ", int(a[i].y/2 + 0.5)); } return 0;}

例题：P1919 【模板】FFT快速傅里叶变换

///优化版#include using namespace std;const int N = 5e6 + 7;const double PI = acos(-1);int n, m;int limit = 1;int res, ans[N];int l;int r[N];struct Complex{ double x, y; Complex (double x = 0, double y = 0) : x(x), y(y) {}} a[N], b[N];Complex operator * (Complex J, Complex Q){ //模长相乘，幅度相加 return Complex(J.x * Q.x - J.y * Q.y, J.x * Q.y + J.y * Q.x);}Complex operator + (Complex J, Complex Q){ return Complex(J.x + Q.x, J.y + Q.y);}Complex operator - (Complex J, Complex Q){ return Complex(J.x - Q.x, J.y - Q.y);}void FFT(Complex *A, int type){ for(int i = 0; i < limit; i++) { if(i < r[i]) swap(A[i], A[r[i]]); ///保证不在换过去 } ///从底层开始合并： for(int mid = 1; mid < limit; mid = mid * 2) { ///待合并区间长度的一半，最开始是长度为1的合并，类似倍增的思想 Complex wn(cos(PI / mid), type * sin(PI / mid)); //单位根 for(int len = mid * 2, pos = 0; pos < limit; pos += len) { ///len是区间的长度，pos是当前的位置 Complex w(1, 0); for(int k = 0; k < mid; k++, w = w * wn) { ///只扫左半部分，蝴蝶变换得到有半部分。 Complex x = A[pos + k]; //左半部分 Complex y = w * A[pos + mid + k]; //有半部分 A[pos + k] = x + y; A[pos + mid + k] = x - y; } } } if(type == 1) return ; for(int i = 0; i <= limit; i++) { ans[i] += (int)(A[i].y / limit/2 + 0.5); if(ans[i] >= 10) ///考虑进位问题。 { ans[i + 1] = ans[i] / 10; ans[i] %= 10; limit += (i == limit); ///进位 } ///除以我们推出的N。 }}string str, s;int main(){ cin >> str >> s; n = str.size(); m = s.size(); int cnt, tot; cnt = tot = 0; for(int i = n - 1; i >= 0; i--) a[cnt++].x = str[i] - '0'; for(int i = m - 1; i >= 0; i--) a[tot++].y = s[i] - '0'; while(limit < n + m) limit <<= 1, l++; for(int i = 0; i <= limit; i++) { r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); } FFT(a, 1); ///FFT(b, 1); for(int i = 0; i <= limit; ++i) { a[i] = a[i] * a[i]; } FFT(a, -1); while(!ans[limit] && limit >= 1) limit--; ///去除前导0 limit++; while(--limit >= 0) { cout << ans[limit]; } return 0;}

正常版。

//正常的#include using namespace std;const int N = 5e6 + 7;const double PI = acos(-1);int n, m;int limit = 1;int res, ans[N];int l;int r[N];struct Complex{ double x, y; Complex (double x = 0, double y = 0) : x(x), y(y) {}} a[N], b[N];Complex operator * (Complex J, Complex Q){ //模长相乘，幅度相加 return Complex(J.x * Q.x - J.y * Q.y, J.x * Q.y + J.y * Q.x);}Complex operator + (Complex J, Complex Q){ return Complex(J.x + Q.x, J.y + Q.y);}Complex operator - (Complex J, Complex Q){ return Complex(J.x - Q.x, J.y - Q.y);}void FFT(Complex *A, int type){ for(int i = 0; i < limit; i++) { if(i < r[i]) swap(A[i], A[r[i]]); ///保证不在换过去 } ///从底层开始合并： for(int mid = 1; mid < limit; mid = mid * 2) { ///待合并区间长度的一半，最开始是长度为1的合并，类似倍增的思想 Complex wn(cos(PI / mid), type * sin(PI / mid)); //单位根 for(int len = mid * 2, pos = 0; pos < limit; pos += len) { ///len是区间的长度，pos是当前的位置 Complex w(1, 0); for(int k = 0; k < mid; k++, w = w * wn) { ///只扫左半部分，蝴蝶变换得到有半部分。 Complex x = A[pos + k]; //左半部分 Complex y = w * A[pos + mid + k]; //有半部分 A[pos + k] = x + y; A[pos + mid + k] = x - y; } } } if(type == 1) return ; for(int i = 0; i <= limit; i++) { ans[i] += (int)(A[i].x / limit + 0.5); if(ans[i] >= 10) ///考虑进位问题。 { ans[i + 1] = ans[i] / 10; ans[i] %= 10; limit += (i == limit); ///进位 } ///除以我们推出的N。 }}string str, s;int main(){ cin >> str >> s; n = str.size(); m = s.size(); int cnt, tot; cnt = tot = 0; for(int i = n - 1; i >= 0; i--) a[cnt++].x = str[i] - '0'; for(int i = m - 1; i >= 0; i--) b[tot++].x = s[i] - '0'; while(limit < n + m) limit <<= 1, l++; for(int i = 0; i <= limit; i++) { r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); } FFT(a, 1); FFT(b, 1); for(int i = 0; i <= limit; ++i) { a[i] = a[i] * b[i]; } FFT(a, -1); while(!ans[limit] && limit >= 1) limit--; ///去除前导0 limit++; while(--limit>=0) { cout << ans[limit]; } return 0;}

P3338 [ZJOI2014]力

先放上，防止自己咕咕咕

题意：

n个数

快速数论变换 (NTT)

#include using namespace std;#define ll long longconst int N = 5e6 + 7;const double PI = acos(-1);const int p = 998244353, G = 3, Gi = 332748118;int n, m;int limit = 1;int res, ans[N];int l;int r[N];ll a[N], b[N];ll qpow(ll a, ll b){ ll ans = 1; while(b) { if(b & 1) ans = ans * a % p; a = a * a % p; b >>= 1; } return ans;}void NTT(ll *A, int type){ for(int i = 0; i < limit; i++) { if(i < r[i]) swap(A[i], A[r[i]]); ///保证不在换过去 } ///从底层开始合并： for(int mid = 1; mid < limit; mid = mid * 2) { ///待合并区间长度的一半，最开始是长度为1的合并，类似倍增的思想 ll wn = qpow(G, (p - 1) / (mid * 2)); if(type==-1) wn = qpow(wn, p - 2); for(int len = mid * 2, pos = 0; pos < limit; pos += len) { ///len是区间的长度，pos是当前的位置 ll w = 1; for(int k = 0; k < mid; k++, w = w * wn % p) { ///只扫左半部分，蝴蝶变换得到有半部分。 int x = A[pos + k]; //左半部分 int y = w * A[pos + mid + k] % p; //有半部分 A[pos + k] = (x + y) % p; A[pos + mid + k] = (x - y + p) % p; } } } if(type == 1) return ; ll inlimit = qpow(limit, p - 2); for(int i = 0; i < limit; i++) { A[i] = (A[i] * inlimit) % p; ///除以我们推出的N。 }}int main(){ cin >> n >> m; for(int i = 0; i <= n; i++) cin >> a[i], a[i] = (a[i] + p) % p; for(int i = 0; i <= m; i++) cin >> b[i], b[i] = (b[i] + p) % p; while(limit <= n + m) limit <<= 1, l++; for(int i = 0; i <= limit; i++) { r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); } NTT(a, 1); NTT(b, 1); for(int i = 0; i <= limit; ++i) { a[i] = a[i] * b[i] % p; } NTT(a, -1); for(int i = 0; i <= n + m; i++) cout << a[i] << " "; return 0;}

版权声明：本文内容由网络用户投稿，版权归原作者所有，本站不拥有其著作权，亦不承担相应法律责任。如果您发现本站中有涉嫌抄袭或描述失实的内容，请联系我们jiasou666@gmail.com 处理，核实后本网站将在24小时内删除侵权内容。