高精度乘法求调
查看原帖
高精度乘法求调
482610
Mortidesperatslav楼主2024/1/21 12:26

rt,我用 P1720 做 FFT+矩阵快速幂斐波那契数列板题。

#include<bits/stdc++.h>
using namespace std;
struct mycomplex {
	double x,y;
	mycomplex(double tmpx=0,double tmpy=0) {
		this->x=tmpx,this->y=tmpy;
	}
	mycomplex operator+(mycomplex b) {
		return {this->x+b.x,this->y+b.y};
	}
	mycomplex operator-(mycomplex b) {
		return {this->x-b.x,this->y-b.y};
	}
	mycomplex operator*(mycomplex b) {
		return {this->x*b.x-this->y*b.y,this->x*b.y+this->y*b.x};
	}
}a[50005],b[50005];
int n,m,l,r[50005],ans[50005];
int lim=1;
double pi=3.1415926535897932;
void FFT(mycomplex *aa,int tp) {
	for(int i=0; i<lim; i++)if(i<r[i])swap(aa[i],aa[r[i]]);
	for(int mid=1; mid<lim; mid<<=1) {
		mycomplex omega(cos(pi/mid),tp*sin(pi/mid));
		for(int rr=mid<<1,j=0; j<lim; j+=rr) {
			mycomplex omega2(1,0);
			for(int k=0; k<mid; k++,omega2=omega2*omega) {
				mycomplex xx=aa[j+k],yy=omega2*aa[j+k+mid];
				//	cout<<xx.x<<" "<<xx.y<<" "<<yy.x<<" "<<yy.y<<"\n";
				aa[j+k]=xx+yy;
				aa[j+k+mid]=xx-yy;
				//	cout<<aa[j+k].x<<" "<<aa[j+k].y<<"\n";
			}
		}
	}
}
class bigint_256 {
	private:
		int val[30005];
		int size;
	public:
		bigint_256() {
			memset(this->val, 0, sizeof(this->val));
		}
		friend std::istream& operator >>(std::istream &is, bigint_256 &x) { //input
			std::string tmp;
			is >> tmp;
			int len = tmp.length();
			for (register int i = 0; i < len; i++)x.val[i] = tmp[len - i - 1] - '0';
			x.size = len;
			return is;
		}
		friend std::ostream& operator <<(std::ostream &os, bigint_256 x) { //output
			for (register int i = x.size - 1; i >= 0; i--)os << x.val[i];
			return os;
		}
		void operator =(bigint_256 b) {
			for (register int i = 0; i < b.size; i++)this->val[i] = b.val[i];
			this->size = b.size;
		}
		void operator =(string b) {
			for (register int i = 0; i < b.size(); i++)this->val[i] = b[i]-'0';
			this->size = b.size();
		}
		bool operator ==(bigint_256 b) {
			if(this->size != b.size)return 0;
			for (register int i = 0; i < b.size; i++)if(this->val[i] != b.val[i])return 0;
			return 1;
		}
		bigint_256 operator +(bigint_256 b) {
			bigint_256 tmp=b;
			short jw=0;
			for(register int i=0; i<=max(this->size,tmp.size)+1; i++) {
				tmp.val[i]=this->val[i]+b.val[i]+jw;
				jw=tmp.val[i]/10;
				tmp.val[i]%=10;
			}
			tmp.size=max(this->size,tmp.size)+1;
			while(tmp.val[tmp.size-1]==0&&tmp.size>1)tmp.size--;
			return tmp;
		}
		void operator +=(bigint_256 b) {
			short jw = 0;
			for (register int i = 0; i <= max(this->size, b.size) + 1; i++) {
				this->val[i] += (b.val[i] + jw);
				jw = this->val[i] / 10;
				this->val[i] %= 10;
			}
			this->size = max(this->size, b.size) + 1;
			while (this->val[this->size - 1] == 0 && this->size > 1)this->size--;
		}
		bigint_256 operator -(bigint_256 b) {
			bigint_256 tmp = b;
			short jw = 0;
			for (register int i = 0; i <= max(this->size, tmp.size) + 1; i++) {
				tmp.val[i] = this->val[i] - b.val[i] + jw;
				if (tmp.val[i] < 0)jw = -1;
				else jw = 0;
				tmp.val[i] = (tmp.val[i] + 1000) % 10;
			}
			tmp.size = max(this->size, tmp.size) + 1;
			while (tmp.val[tmp.size - 1] == 0 && tmp.size > 1)tmp.size--;
			return tmp;
		}
		void operator -=(bigint_256 b) {
			short jw = 0;
			for (register int i = 0; i <= max(this->size, b.size) + 1; i++) {
				this->val[i] -= (b.val[i] + jw);
				jw = this->val[i] / 10;
				this->val[i] = (this->val[i] + 1000) % 10;
			}
			this->size = max(this->size, b.size) + 1;
			while (this->val[this->size - 1] == 0 && this->size > 1)this->size--;
		}
		bigint_256 operator *(bigint_256 bb) {
			bigint_256 ans;
			n=this->size-1;
			m=bb.size-1;
			l=0,lim=1;
			for(int i=0;i<=n+m;i++)a[i].x=b[i].x=0;
			for(int i=0; i<=n; i++)a[i].x=this->val[i];
			for(int i=0; i<=m; i++)b[i].x=bb.val[i];
			while(lim<=n+m)lim<<=1,l++;
			memset(r,0,sizeof(r));
			for(int i=0; i<lim; i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
			FFT(a,1);
			FFT(b,1);
			for(int i=0; i<=lim; i++)a[i]=a[i]*b[i];
			FFT(a,-1);
			for(int i=0; i<=n+m; i++)ans.val[i]=(int)(a[i].x/lim+0.5);
			for(int i=0; i<=n+m+1; i++)
				if(ans.val[i]>=10) {
					ans.val[i+1]+=ans.val[i]/10;
					ans.val[i]%=10;
				}
			ans.size=n+m+2;
			while(ans.val[ans.size-1]==0&&ans.size>1)ans.size--;
			return ans;
		}
};
struct juzhen {
	public:
		bigint_256 a[2][2];
};
juzhen p;
juzhen juzhenchengfa(juzhen x,juzhen y) {
	juzhen res;
	string zr="0";
	res.a[0][0]=zr;
	res.a[1][0]=zr;
	res.a[0][1]=zr;
	res.a[1][1]=zr;
	for(register int i=0; i<2; i++)
		for(register int j=0; j<2; j++)
			for(register int k=0; k<2; k++)
				res.a[i][j]=(res.a[i][j]+(x.a[i][k]*y.a[k][j]));
	return res;
}
juzhen juzhenkuaisumi(juzhen x,int y) {
	juzhen res;
	bool fl=0;
	while(y>0) {
		if(y&1) {
			if(fl==0) {
				res=x;
				fl=1;
			} else res=juzhenchengfa(res,x);
		}
		x=juzhenchengfa(x,x);
		y>>=1;
	}
	return res;
}
int t;
int main() {
	ios::sync_with_stdio(false);
	cin.tie(0);
	cout.tie(0);
//	cin>>t;
//	while(t--) {
		int kk;
		cin>>kk;
		if(kk==0)cout<<"0.00\n";
		else if(kk==1)cout<<"1.00\n";
		else if(kk==2)cout<<"1.00\n";
		else {
		//	if((kk&1)==0)cout<<"-";
			p.a[0][0]="1";
			p.a[0][1]="1";
			p.a[1][0]="1";
			p.a[1][1]="0";
			juzhen q=juzhenkuaisumi(p,kk-2);
			bigint_256 ans=q.a[0][1]+q.a[0][0];
			cout<<ans<<".00"<<"\n";
		}
//	}
}

但是输入 1616 时输出 135-63.00。

不知道什么原因,求调。

2024/1/21 12:26
加载中...