FFT板子RE求调
  • 板块灌水区
  • 楼主godmoo
  • 当前回复5
  • 已保存回复5
  • 发布时间2024/7/25 19:56
  • 上次更新2024/7/25 20:46:47
查看原帖
FFT板子RE求调
1109707
godmoo楼主2024/7/25 19:56
#include<bits/stdc++.h>
#define pb push_back
#define is insert
#define fi first
#define se second
#define bg begin
#define INF INT_MAX
#define Mod(a,m) (((a)%(m)+(m))%(m))
#define mem(a) memset(a,0,sizeof a)
#define cpy(a,b) memcpy(a,b,sizeof b)
//#pragma comment(linker,"/STACK:1024000000,1024000000")
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
const int N=1<<21;
const double PI=acos(-1);
int n,m;
vector<int> a,b;
class Complex{
public:
	double a,b;
	Complex(double _a=0.0,double _b=0.0):a(_a),b(_b){}
	Complex operator +(const Complex _) const{return Complex(a+_.a,b+_.b);}
	Complex operator -(const Complex _) const{return Complex(a-_.a,b-_.b);}
	Complex operator *(const Complex _) const{return Complex(a*_.a-b*_.b,a*_.b+b*_.a);}
};
typedef Complex cpd;
class Poly{
private:
	int st;
	vector<cpd> v;
public:
	Poly(int l,int _st){st=_st;for(int i=1;i<=l;i++)putchar('f'),v.pb(Complex(0,0));}
	Poly(int* a,int l){v.resize(l);st=-1;for(int i=0;i<l;i++)v[i]=Complex(a[i],0);}
	Poly(vector<int> a){v.resize(a.size());st=-1;for(int i=0;i<v.size();i++)v[i]=Complex(a[i],0);}
	Poly(const Poly &_){st=_.st;v.clear();for(auto c:_.v)v.pb(c);} 
	int size(){return v.size();}
	void resize(int l){v.resize(l);}
	int suit(){return (1<<__lg(v.size()))==v.size()?v.size():(1<<(__lg(v.size())+1));}
	double& operator [](int _){assert(_<v.size());if(~st)idft(size());return v[_].a;}
	vector<int> formula(){cout<<"before formula!\n";if(~st)idft(size());cout<<"after judge!\n";vector<int> res;cout<<"create vector!\n";for(auto c:v)res.pb((int)c.a);cout<<"perfect!\n";return res;}
	Poly operator +(const Poly _) const{
		Poly A(*this),B(_);
		int l=max(A.suit(),B.suit());
		A.resize(l),B.resize(l);
		Poly res(l,A.st);
		for(int i=0;i<l;i++) res.v[i]=A.v[i]+B.v[i];
		return res;
	}
	Poly operator -(const Poly _) const{
		Poly A(*this),B(_);
		int l=max(A.suit(),B.suit());
		A.resize(l),B.resize(l);
		Poly res(l,A.st);
		for(int i=0;i<l;i++) res.v[i]=A.v[i]-B.v[i];
		return res;
	}
private:
	int rev[N];
	void bitrev(int l){
		for(int i=0;i<l;i++){
			rev[i]=rev[i>>1]>>1;
			if(i&1) rev[i]|=l>>1;
		}
		for(int i=0;i<l;i++) if(i<rev[i]) swap(v[i],v[rev[i]]);
	}
	void fft(int l,int op){
		assert(st!=op);
		st=op;
		cout<<"before bitrev!\n";
		bitrev(l);
		cout<<"before merge!\n";
		for(int i=2;i<=l;i<<=1){
			cpd omega(cos(2*PI/i),sin(op*2*PI/i));
			for(int j=0;j<l;j+=i){
				cpd x(1,0);
				for(int k=j;k<=j+(i>>1);k++){
					cpd g=v[k],xh=x*v[k+(i>>1)];
					v[k]=g+xh;
					v[k+(i>>1)]=g-xh;
					x=x*omega;
				}
			}
		}
	}
	void dft(int l){fft(l,1);}
	void idft(int l){fft(l,-1);for(int i=0;i<l;i++)v[i].a/=l;}
public:
	Poly operator *(const Poly _) const{
		Poly A(*this),B(_);
		if(~A.st) A.idft(A.size());
		if(~B.st) B.idft(B.size());
		cout<<"prepare!\n";
		int s=A.size()+B.size()-1;
		int l=(1<<__lg(s))==s?s:(1<<(__lg(s)+1));
		A.resize(l),B.resize(l);
		cout<<"before dft!\n";
		A.dft(l),B.dft(l);
		cout<<"end dft!\n" ;
		cout<<l<<endl;
		Poly res(l,1);
		cout<<"bwrt";
		for(int i=0;i<l;i++) res.v[i]=A.v[i]*B.v[i];
		cout<<"before idft!\n";
		res.idft(l);
		cout<<"end idft!\n";
		return res;
	}
};
int main() {
    // 初始化两个多项式 1x^2 + 2x + 3 和 4x^2 + 5x + 6
    vector<int> coeffs1 = {1, 2, 3}; // 对应1x^2 + 2x + 3
    vector<int> coeffs2 = {4, 5, 6}; // 对应4x^2 + 5x + 6

    Poly poly1(coeffs1);
    Poly poly2(coeffs2);

    // 多项式加法
    Poly poly_add = poly1 + poly2;
    cout << "Polynomial Addition: ";
    for (auto c : poly_add.formula()) {
        cout << c << " ";
    }
    cout << endl;

    // 多项式减法
    Poly poly_sub = poly1 - poly2;
    cout << "Polynomial Subtraction: ";
    for (auto c : poly_sub.formula()) {
        cout << c << " ";
    }
    cout << endl;

    // 多项式乘法
    Poly poly_mul = poly1 * poly2;
    cout << "Polynomial Multiplication: ";
    for (auto c : poly_mul.formula()) {
        cout << c << " ";
    }
    cout << endl;

    return 0;
}

2024/7/25 19:56
加载中...