FFT样例过 0分RE求助(带注释)
查看原帖
FFT样例过 0分RE求助(带注释)
399475
_XHY20180718_楼主2022/9/16 13:52

RE,

#include<cstdio>
#include<cmath>
#include<algorithm>
typedef long long int ll;
using namespace std;
const ll MOD=1e8;
const int MAXN=2^18;

inline void readp(ll *num)//高精度从高位到低位读入
{
	ll tmp[MAXN];//从高位到低位暂存
	register int i=0;//十进制位数
	register int m=0,n=0;//八位压位十进值位数
	register ll s=0,t=0;//暂存压位位值
	register char ch=getchar();
	while(ch<'0'||ch>'9') ch=getchar();	
	while(ch>='0'&&ch<='9') 
	{
		i++;
		s=(s<<3)+(s<<1)+(ch-'0');
		if(i%8==0)
		{
			tmp[++m]=s;
			s=0;			
		}
		ch=getchar();	
	}
	
	int k=i%8;//恢复运算逆序存储
	if(k==0) 
	{
		num[0]=m;
		for(; m>=1; m--)
		{
			num[++n]=tmp[m];	
		}
		return;
	}
	//当前现有k位加上前8位的8-k组成后8位
	ll sl=1;
	while(k>=1)//最后压位数
	{
		sl*=10;
		k--;
	}
	ll sr=MOD/sl;//sl=10^k,sr=10^(8-k)
	tmp[0]=0;
	for(; m>=0; m--)
	{
		t=tmp[m]%sr;//前8位的低8-k位
		num[++n]=s+t*sl;//后8位组成:前8位的低8-k位和当前现有k位
		s=tmp[m]/sr;//前8位的高k位
	}	
	num[0]=n;
	return;
} 

inline void writei(ll s)//高精度从高位到低位写出 
{
	if(s>=10) writei(s/10);
	putchar((s%10)+'0');	
}

//int整数数组高精度
inline void writep(ll num[])//高精度从高位到低位写出 
{
	register int i=0;
	for(i=num[0]; i>=1; i--) writei(num[i]);	
}

struct node//复数结构体
{
    double x, y;
    node(double xx=0, double yy=0)//复数初始化实部x与虚部y为0
	{
      x=xx, y=yy;
    }
}A[MAXN], B[MAXN], C[MAXN];

//重定义复数运算
node operator *(node z1, node z2)
{
  return node(z1.x*z2.x-z1.y*z2.y, z1.x*z2.y+z1.y*z2.x);
}

node operator +(node z1, node z2)
{
  return node(z1.x+z2.x, z1.y+z2.y);
}

node operator -(node z1, node z2)
{
  return node(z1.x-z2.x, z1.y-z2.y);
}

const double Pi=acos(-1.0);
int FL=1, L=0;//FFT长度和长度二进制2^L=FL
int R[MAXN];//蝴蝶递推反序

inline void fft(node *zs, double fl)
{
	for(register int i=0; i<FL; i++)//蝴蝶递推反序
	{
   	if(i<R[i]) swap(zs[i],zs[R[i]]);//前面的if保证只换一次
	}  
	for(register int i=1; i<FL; i<<=1)//迭代蝴蝶递推
	{
  	node T(cos(Pi/i), fl*sin(Pi/i));//单位根
    for(register int j=0; j<FL; j+=(i<<1))
		{
      node t(1, 0);
      for(register int k=0; k<i; k++, t=t*T)
			{
      	node Nx=zs[j+k];
      	node Ny=t*zs[i+j+k];
       	zs[j+k]=Nx+Ny;
       	zs[i+j+k]=Nx-Ny;
   		}
    }
	}
}

ll a[MAXN];//乘数 
ll b[MAXN];//乘数
ll c[MAXN];//积 

inline void mul_hc(const ll s[], const ll t[], ll *r)
{
	//乘法高精度需要从低位到高位运算

	//读入的数的每一位看成多项式的一项,保存在复数的实部 
	register int N=0,M=0;//多项式次数
  for(register int i=1; i<=s[0]; i++) A[i-1].x=(double)(s[i]), N++;
	for(register int i=1; i<=t[0]; i++) B[i-1].x=(double)(t[i]), M++;
	//register int FL=1,L=0;//FFT长度和长度二进制2^L=FL
	while(FL<N+M) FL<<=1, L++;
	//FFT二进制反序:蝴蝶递推定理:例如:FL=8 L=3
	//原来的序号:0、1、2、3、4、5、6、7
	//现在的序号:0、4、2、6、1、5、3、7
	for(register int i=0; i<FL; i++) 
	{
		R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
	}
	fft(A, 1);
	fft(B, 1);
  for(register int i=0; i<=FL; i++) C[i]=A[i]*B[i];//时域卷积等价频域乘积
  fft(C, -1);
	for(register int i=0; i<=FL; i++) 
	{
  	r[i+1]+=(ll)(C[i].x/FL+0.5);
   	if(r[i+1]>=MOD)
		{
		 	r[i+2]+=r[i+1]/MOD;
			r[i+1]%=MOD;
			FL+=(i==FL);
		}
  }
	while(r[FL+1]==0&&FL>=0) FL--;
	r[0]=FL+1;
}

int main()
{
	readp(a);
	readp(b);
	mul_hc(a,b,c); 
	writep(c);
	putchar('\n'); 
	return 0;
}

2022/9/16 13:52
加载中...