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