4fd800b3a8 2011-02-23        kinaba: 
4fd800b3a8 2011-02-23        kinaba: //-------------------------------------------------------------
4fd800b3a8 2011-02-23        kinaba: // Fast Discrete Fourier Transform
4fd800b3a8 2011-02-23        kinaba: //   O(n log n)
4fd800b3a8 2011-02-23        kinaba: //   n must be a power of 2.
4fd800b3a8 2011-02-23        kinaba: //
4fd800b3a8 2011-02-23        kinaba: // Verified by
4fd800b3a8 2011-02-23        kinaba: //   - SRM 436 LV3
4fd800b3a8 2011-02-23        kinaba: //-------------------------------------------------------------
4fd800b3a8 2011-02-23        kinaba: 
4fd800b3a8 2011-02-23        kinaba: CMP tmp[65536*4];
4fd800b3a8 2011-02-23        kinaba: 
4fd800b3a8 2011-02-23        kinaba: template<int F>
4fd800b3a8 2011-02-23        kinaba: void fft_impl( CMP a[], int n, int stride = 1 )
4fd800b3a8 2011-02-23        kinaba: {
4fd800b3a8 2011-02-23        kinaba: 	if( n > 1 )
4fd800b3a8 2011-02-23        kinaba: 	{
4fd800b3a8 2011-02-23        kinaba: 		CMP *ev=a, *od=a+stride;
4fd800b3a8 2011-02-23        kinaba: 		int s2=stride*2, n2=n/2;
4fd800b3a8 2011-02-23        kinaba: 
4fd800b3a8 2011-02-23        kinaba: 		fft_impl<F>( ev, n2, s2 );
4fd800b3a8 2011-02-23        kinaba: 		fft_impl<F>( od, n2, s2 );
4fd800b3a8 2011-02-23        kinaba: 
4fd800b3a8 2011-02-23        kinaba: 		for(int i=0; i<n; ++i) tmp[i] = ev[s2*(i%n2)] + od[s2*(i%n2)]*polar(1.0, F*2*M_PI*i/n);
4fd800b3a8 2011-02-23        kinaba: 		for(int i=0; i<n; ++i) a[stride*i] = tmp[i];
4fd800b3a8 2011-02-23        kinaba: 	}
4fd800b3a8 2011-02-23        kinaba: }
4fd800b3a8 2011-02-23        kinaba: 
4fd800b3a8 2011-02-23        kinaba: void fft( vector<CMP>& a )
4fd800b3a8 2011-02-23        kinaba: {
4fd800b3a8 2011-02-23        kinaba: 	fft_impl<+1>(&a[0], a.size());
4fd800b3a8 2011-02-23        kinaba: }
4fd800b3a8 2011-02-23        kinaba: 
4fd800b3a8 2011-02-23        kinaba: void ifft( vector<CMP>& a )
4fd800b3a8 2011-02-23        kinaba: {
4fd800b3a8 2011-02-23        kinaba: 	fft_impl<-1>(&a[0], a.size());
4fd800b3a8 2011-02-23        kinaba: 	for(int i=0; i<a.size(); ++i)
4fd800b3a8 2011-02-23        kinaba: 		a[i] /= a.size();
4fd800b3a8 2011-02-23        kinaba: }