FFT(快速博立叶变换)

xiaoxiao2021-02-28  20

证明请参考算法概论的快速博立叶。

FFT递归:

#include <bits/stdc++.h> using namespace std; const int maxn = 5005; const double PI = acos ( -1.0 ); struct node { double real, image; node ( double a = 0, double b = 0 ) : real ( a ), image ( b ) { } node operator + ( const node &b ) { return node ( real+b.real, image+b.image ); } node operator - ( const node &b ) { return node ( real-b.real, image-b.image ); } node operator * ( const node &b ) { return node ( real*b.real-image*b.image, real*b.image+image*b.real ); } }; node a[maxn], b[maxn]; int re[maxn]; node * FFT ( node x[maxn], int n, int dft ) { node * ans = ( node * ) malloc ( sizeof ( node )*n ); if ( n == 1 ) { ans[0] = x[0]; return ans; } node lf[maxn], rt[maxn]; for ( int i = 0; i < n; i += 2 ) lf[i/2] = x[i]; for ( int i = 1; i < n; i += 2 ) rt[i/2] = x[i]; node * L = FFT ( lf, n/2, dft ); node * R = FFT ( rt, n/2, dft ); node nw = node ( cos ( dft*2*PI/n ), sin ( dft*2*PI/n ) ); node w = node ( 1, 0 ); for ( int i = 0; i < n/2; i ++ ) { ans[i] = L[i]+w*R[i]; ans[i+n/2] = L[i]-w*R[i]; w = w*nw; } return ans; } void solve ( ) { int n; scanf ( "%d", &n ); for ( int i = 0; i < n; i ++ ) { scanf ( "%lf", &a[i].real ); a[i].image = 0; } for ( int i = 0; i < n; i ++ ) { scanf ( "%lf", &b[i].real ); b[i].image = 0; } int p = 1; while ( p < 2*n ) p <<= 1; for ( int i = n; i < p; i ++ ) a[i].real = a[i].image = b[i].real = b[i].image = 0; node * pa = FFT ( a, p, 1 ); node * pb = FFT ( b, p, 1 ); for ( int i = 0; i < p; i ++ ) a[i] = pa[i]*pb[i]; pa = FFT ( a, p, -1 ); for ( int i = 0; i < p; i ++ ) re[i] = pa[i].real/p+0.5; for ( int i = 0; i < p; i ++ ) printf ( "%d ", re[i] ); printf ( "\n" ); } int main ( ) { solve ( ); return 0; }

FFT迭代代码(HDU1402):

#include <bits/stdc++.h> using namespace std; const int maxn = 500005; const double PI = acos ( -1.0 ); struct node { double real, image; node ( double a = 0, double b = 0 ) : real ( a ), image ( b ) { } node operator + ( const node &b ) { return node ( real+b.real, image+b.image ); } node operator - ( const node &b ) { return node ( real-b.real, image-b.image ); } node operator * ( const node &b ) { return node ( real*b.real-image*b.image, real*b.image+image*b.real ); } }; node a[maxn], b[maxn], W[maxn]; int re[maxn]; char x[maxn], y[maxn]; void change ( node a[], int n ) { int p = 1; while ( ( 1 << p ) < n ) p ++; for ( int i = 0; i < n; i ++ ) { int j = 0, k = i, t = p; while ( t -- ) { j <<= 1; j |= k&1; k >>= 1; } if ( j > i ) swap ( a[i], a[j] ); } } void FFT ( node a[], int n, int dft ) { change ( a, n ); for ( int i = 1; i < n; i *= 2 ) { node nw ( cos ( dft*PI*2/2/i ), sin ( dft*PI*2/2/i ) ); for ( int j = 0; j < n; j += 2*i ) { node w ( 1, 0 ); for ( int k = 0; k < i; k ++ ) { node tmp = a[j+k]; node res = w*a[j+k+i]; a[j+k] = tmp+res; a[j+k+i] = tmp-res; w = w*nw; } } } if ( dft == -1 ) { for ( int i = 0; i < n; i ++ ) a[i].real /= n; } } void solve ( ) { int lenx, leny; while ( ~ scanf ( "%s%s", x, y ) ) { lenx = strlen ( x ); for ( int i = 0; i < lenx; i ++ ) { a[i].real = x[lenx-1-i]-'0'; a[i].image = 0; } leny = strlen ( y ); for ( int i = 0; i < leny; i ++ ) { b[i].real = y[leny-1-i]-'0'; b[i].image = 0; } int p = 1; while ( p < 2*lenx || p < 2*leny ) p <<= 1; for ( int i = lenx; i < p; i ++ ) a[i].real = a[i].image = 0; for ( int i = leny; i < p; i ++ ) b[i].real = b[i].image = 0; FFT ( a, p, 1 ); FFT ( b, p, 1 ); for ( int i = 0; i < p; i ++ ) a[i] = a[i]*b[i]; FFT ( a, p, -1 ); for ( int i = 0; i < p; i ++ ) re[i] = a[i].real+0.5; for ( int i = 0; i < p; i ++ ) { re[i+1] += re[i]/10; re[i] %= 10; } int h = 0; for ( int i = p-1; i >= 0; i -- ) if ( re[i] ) { h = i; break ; } for ( int i = h; i >= 0; i -- ) printf ( "%d", re[i] ); printf ( "\n" ); } } int main ( ) { solve ( ); return 0; }

转载请注明原文地址: https://www.6miu.com/read-850071.html

最新回复(0)