FFT를 위한 Poly 클래스 (재귀) || Polynomial class for FFT (recursive)

polyfft.cpp

#include<iostream>
#include<complex>
#include<vector>
#include<cmath>
#include<initializer_list>
using namespace std;

#define FLAG_DOROUND 1

const double Pi=acos(-1);

template<class T>
struct cpxProxy
{
	using cpx=complex<T>;
	cpx& orig;
	operator T(){return orig.real();}
	cpx& raw(){return orig;}
	cpxProxy<T>& operator=(T&& n)
	{
		orig=n;return *this;
	}
};
template<class T>
istream& operator>>(istream& is,cpxProxy<T>& cp)
{
	T t;cin>>t;cp=t;return is;
}

template<class T>
struct Poly
{
	using cpx=complex<T>;
	using proxy=cpxProxy<T>;
	vector<cpx> C;size_t S;
	Poly(size_t s):C(s,0),S(s){}
	Poly(initializer_list<cpx> IL):C(IL),S(size(IL)){}
	proxy operator[](size_t&& s)
	{
		return proxy{C[s]};
	}
	void FFT(cpx w)
	{
		if(S==1)return;
		Poly odd(S>>1),even(S>>1);
		for(size_t i=0;i<S;i++)
		(i&1?odd:even).C[i>>1]=C[i];
		odd.FFT(w*w);
		even.FFT(w*w);
		cpx wp(1,0);
		for(size_t i=0;i<S>>1;i++)
		{
			C[i]=even.C[i]+odd.C[i]*wp;
			C[i+(S>>1)]=even.C[i]-odd.C[i]*wp;
			wp*=w;
		}
	}
	Poly& operator*=(Poly P)
	{
		(*this)=(*this)*P;
		return (*this);
	}
	void resize(size_t s){C.resize(s);S=s;}
};
template<class T>
Poly<T> operator*(Poly<T> A, Poly<T> B)
{
	size_t n=1,resz=A.S+B.S;
	while(n<resz)n<<=1;
	A.resize(n);
	B.resize(n);
	Poly<T> C(n);
	complex<T> w(cos(2*Pi/n),sin(2*Pi/n));
	A.FFT(w); B.FFT(w);
	for(size_t i=0;i<n;i++)
	C.C[i]=A.C[i]*B.C[i];
	C.FFT(complex<T>(1,0)/w);
	for(size_t i=0;i<n;i++)
	{
		C.C[i]/=complex<T>(n,0);
		#ifdef FLAG_DOROUND
		C.C[i]=
		complex<T>(round(C.C[i].real()),round(C.C[i].imag()));
		#endif
	}
	C.resize(resz);
	return C;
}