Eps开到几,平生我自知。为什么Eps要开到50000啊?!
查看原帖
Eps开到几,平生我自知。为什么Eps要开到50000啊?!
658786
STUDENT00楼主2023/12/21 19:29

RT,FFT,double判精度这么菜!Eps开到50000才过?!

#include<bits/stdc++.h>
#define complex complex<double>
using namespace std;
const int N=1200005;
const double PI=acos(-1);
int n,m;
int A[N],B[N];
char Sa[N],Sb[N];
complex a[N],b[N],p[N];
/*void FFT(complex *A,int n,int op){
	if(n==1) return;
	complex A1[n/2],A2[n/2];
	for(int i=0;i<n/2;i++) A1[i]=A[2*i],A2[i]=A[2*i+1];
	FFT(A1,n/2,op),FFT(A2,n/2,op);
	complex w1({cos(2*PI/n),sin(2*PI/n)*op});
	complex wk({1,0});
	for(int i=0;i<n/2;i++){
		A[i]=A1[i]+A2[i]*wk;
		A[i+n/2]=A1[i]-A2[i]*wk;
		wk*=w1;
	}
}*/
void FFT(complex *A,int n,int op){
	if(n==1) return;
	complex A1[n/2],A2[n/2];
	for(int i=0;i<n/2;i++) A1[i]=A[2*i],A2[i]=A[2*i+1];
	FFT(A1,n/2,op),FFT(A2,n/2,op);
	complex w1({cos(2*PI/n),sin(2*PI/n)*op});
	complex wk({1,0});
	for(int i=0;i<n/2;i++){
		A[i]=A1[i]+A2[i]*wk;
		A[i+n/2]=A1[i]-A2[i]*wk;
		wk*=w1;
	}
}
int main(){
	complex two=2;
	scanf("%d%d",&m,&n);
	int tmp=n;
	int k;
	for(k=tmp<<1,n=1;n<=k;n<<=1);
	scanf("%s%s",Sa,Sb);
	for(int i=0;i<m;i++) A[i]=(Sa[i]=='*'?0:Sa[i]-'a'+1);
	for(int i=0;i<tmp;i++) B[i]=(Sb[i]=='*'?0:Sb[i]-'a'+1);
	reverse(A,A+m); 
	//for(int i=0;i<m;i++) cout<<A[i]<<" ";cout<<endl;
	//for(int i=0;i<tmp;i++) cout<<B[i]<<" ";cout<<endl;
	/*for(int i=0;i<m;i++) cout<<A[i];cout<<endl;
	for(int i=0;i<tmp;i++) cout<<B[i];cout<<endl;*/
	for(int i=0;i<n;i++) a[i]=A[i]*A[i]*A[i];
	for(int i=0;i<n;i++) b[i]=B[i];
	FFT(a,n,1),FFT(b,n,1); 
	for(int i=0;i<n;i++) p[i]=a[i]*b[i];
	for(int i=0;i<n;i++) a[i]=A[i];
	for(int i=0;i<n;i++) b[i]=B[i]*B[i]*B[i];
	FFT(a,n,1),FFT(b,n,1); 
	for(int i=0;i<n;i++) p[i]+=a[i]*b[i];
	for(int i=0;i<n;i++) a[i]=A[i]*A[i];
	for(int i=0;i<n;i++) b[i]=B[i]*B[i];
	FFT(a,n,1),FFT(b,n,1); 
	for(int i=0;i<n;i++) p[i]-=two*a[i]*b[i];
	FFT(p,n,-1);
	//for(int i=0;i<n;i++) cout<<p[i].real()<<" ";cout<<endl;
	int cur=0;
	for(int i=0;i<=tmp-m;i++){
		if(fabs(p[i+m-1])<50000/*神秘Eps*/) cur++;
	}
	printf("%d\n",cur);
	for(int i=0;i<=tmp-m;i++){
		if(fabs(p[i+m-1])<50000/*神秘Eps*/) printf("%d ",i+1);
	}
	return 0;
}
2023/12/21 19:29
加载中...