关于FFT
  • 板块学术版
  • 楼主QWQ_HY_DFX
  • 当前回复10
  • 已保存回复11
  • 发布时间2024/1/13 00:15
  • 上次更新2024/1/13 11:48:15
查看原帖
关于FFT
750772
QWQ_HY_DFX楼主2024/1/13 00:15

用到的复数模板:

struct cpx{
	double real,imag;
	inline cpx operator+(cpx c){
		return (cpx){real+c.real,imag+c.imag};
	}
	inline cpx operator-(cpx c){
		return (cpx){real-c.real,imag-c.imag};
	}
	inline cpx operator*(cpx c){
		return (cpx){real*c.real-imag*c.imag,real*c.imag+imag*c.real};
	}
	inline void operator*=(cpx c){
		real = real*c.real-imag*c.imag;
		imag = real*c.imag+imag*c.real;
		return;
	}
}

FFT如果这样写:

inline void FFT(cpx f[],int n,int flag){
  while((1 << bit) < n)
		bit++;
 	n = (1 << bit);
	for(int i = 0;i < n;i++){
    	rev[i] = (rev[i>>1]>>1)|((i&1)<<(bit-1));
    	if(i < rev[i])
			swap(f[i],f[rev[i]]);
    }
	for(int l = 1;l < n;l <<= 1){
		cpx ur = (cpx){cos(pi/(double)l),sin(pi/(double)l)*(double)flag};
		for(int i = 0;i < n;i += (l<<1)){
			cpx omega = (cpx){1,0};
			for(int j = 0;j < l;j++){
				cpx x = f[i+j],y = omega*f[i+j+l];
				f[i+j] = x+y;
				f[i+j+l] = x-y;
				omega = omega*ur;	\\注意这行
			}
		}
	}
	if(flag == -1)
		for(int i = 0;i < n;i++)
			f[i].real /= (double)n;
	return;
}

结果就是对的

如果这样:

inline void FFT(cpx f[],int n,int flag){
  while((1 << bit) < n)
		bit++;
 	n = (1 << bit);
	for(int i = 0;i < n;i++){
    	rev[i] = (rev[i>>1]>>1)|((i&1)<<(bit-1));
    	if(i < rev[i])
			swap(f[i],f[rev[i]]);
    }
	for(int l = 1;l < n;l <<= 1){
		cpx ur = (cpx){cos(pi/(double)l),sin(pi/(double)l)*(double)flag};
		for(int i = 0;i < n;i += (l<<1)){
			cpx omega = (cpx){1,0};
			for(int j = 0;j < l;j++){
				cpx x = f[i+j],y = omega*f[i+j+l];
				f[i+j] = x+y;
				f[i+j+l] = x-y;
				omega *= ur;	//注意这行
			}
		}
	}
	if(flag == -1)
		for(int i = 0;i < n;i++)
			f[i].real /= (double)n;
	return;
}

就是错的
这是精度问题吗?

2024/1/13 00:15
加载中...