LOJ164 高精度除法

调得太难受了,但我并不想多说什么。

压位高精

压18位。最大限度地利用long long__int128__float128,能做\(10^{3000}\)

#include<bits/stdc++.h>
using namespace std;

#define IN inline
#define CO const
typedef long long int64;
typedef __int128 int128;
typedef __float128 float128;

template<class T> IN T read(){
	T x=0,w=1;char c=getchar();
	for(;!isdigit(c);c=getchar())if(c=='-') w=-w;
	for(;isdigit(c);c=getchar()) x=x*10+c-'0';
	return x*w;
}
template<class T> IN T&read(T&x){
	return x=read<T>();
}

CO int64 mod=1e18;
typedef vector<int64> inter;

template<> inter read<inter>(){
	inter ans;
	static char str[3100];
	scanf("%s",str);
	int len=strlen(str);
	for(int i=len-1;i>=0;i-=18){
		int64 sum=0;
		for(int j=max(i-17,0);j<=i;++j) sum=sum*10+str[j]-'0';
		ans.push_back(sum);
	}
	return ans;
}
void write(CO inter&a){
	int n=a.size()-1;
	printf("%lld",a.back());
	for(int i=n-1;i>=0;--i) printf("%018lld",a[i]);
}

bool operator<(CO inter&a,CO inter&b){
	if(a.size()!=b.size()) return a.size()<b.size();
	int n=a.size()-1;
	for(int i=n;i>=0;--i)if(a[i]!=b[i]) return a[i]<b[i];
	return 0;
}
bool operator<=(CO inter&a,CO inter&b){
	if(a.size()!=b.size()) return a.size()<b.size();
	int n=a.size()-1;
	for(int i=n;i>=0;--i)if(a[i]!=b[i]) return a[i]<b[i];
	return 1;
}
bool operator>=(CO inter&a,CO inter&b){
	if(a.size()!=b.size()) return a.size()>b.size();
	int n=a.size()-1;
	for(int i=n;i>=0;--i)if(a[i]!=b[i]) return a[i]>b[i];
	return 1;
}

inter operator+(inter a,CO inter&b){
	int n=a.size()-1,m=b.size()-1;
	if(n<m) n=m,a.resize(n+1);
	for(int i=0;i<=m;++i) a[i]+=b[i];
	for(int i=0;i<n;++i)if(a[i]>=mod) ++a[i+1],a[i]-=mod;
	if(a[n]>=mod) a.push_back(1),a[n]-=mod;
	return a;
}
inter operator-(inter a,CO inter&b){ // a>=b
	int n=a.size()-1,m=b.size()-1;
	for(int i=0;i<=m;++i) a[i]-=b[i];
	for(int i=0;i<=n;++i)if(a[i]<0) --a[i+1],a[i]+=mod;
	for(;n>=1 and !a[n];--n) a.pop_back();
	return a;
}
inter operator*(CO inter&a,CO inter&b){
	if(a==inter{0} or b==inter{0}) return inter{0};
	int n=a.size()-1,m=b.size()-1;
	static int128 tmp[200];
	memset(tmp,0,sizeof(int128)*(n+m+2));
	for(int i=0;i<=n;++i)for(int j=0;j<=m;++j) tmp[i+j]+=(int128)a[i]*b[j];
	for(int i=0;i<=n+m;++i)if(tmp[i]>=mod) tmp[i+1]+=tmp[i]/mod,tmp[i]%=mod;
	inter ans(tmp,tmp+n+m+1);
	if(tmp[n+m+1]) ans.push_back(tmp[n+m+1]);
	return ans;
}
float128 estimate(CO inter&a){
	int n=a.size()-1;
	float128 ans=0;
	for(int i=n;i>=0;--i) ans=ans*mod+a[i];
	return ans;
}
inter operator/(inter a,CO inter&b){
	if(a<b) return inter{0};
	int n=a.size()-1,m=b.size()-1;
	inter ans(n-m+1),rem={0};
	for(int i=n;i>=0;--i){
		rem=rem*inter{0,1}+inter{a[i]};
		if(rem>=b){
			ans[i]=estimate(rem)/estimate(b);
			rem=rem-b*inter{ans[i]};
		}
	}
	while(ans.size()>1 and !ans.back()) ans.pop_back();
	return ans;
}

int main(){
	freopen("B.in","r",stdin),freopen("B.out","w",stdout);
	inter n=read<inter>(),ans={0};
	for(inter i={11};i<=n;i=i*inter{10}+inter{1}) ans=ans+n/i;
	write(ans),puts("");
	return 0;
}

高精度除法

给两个正整数\(a, b\),求\(\lfloor a/b \rfloor\)

对于100%的数据,\(1\leq b\leq a < 10^{100000}\)

题解

倪泽堃《理性愉悦: 高精度数值计算》。

如果实现了高精度实数运算, 那么可以直接运用牛顿迭代法求解大整数\(A\)的倒数。假如考虑方程\(Ax − 1 = 0\),那么迭代式为

\[x=x_0-\frac{Ax_0-1}{A}=\frac{1}{A} \]

很显然,它并没有解决问题。但是如果考虑这样的方程\(\frac{1}{x} − A = 0\),那么迭代式为

\[x=x_0-\frac{\frac{1}{x_0}-A}{-\frac{1}{x_0^2}}=2x_0-Ax_0^2 \]

可以看出,这里只有减法和乘法,因此估计一个较好的初值(可取一个比\(\frac{1}{A}\)稍小的数),就能迭代得到\(\frac{1}{A}\),得到 \(\frac{1}{A}\)后,再乘\(B\),取整后调整,即可得到\(⌊\frac{B}{A}⌋\)

但是,这里用到了不希望出现的高精度实数运算,那么需要寻找一种替代方法。下面提供一种实现的思路。

假如\(A\)\(n\)位, 那么我们希望求得\(A′ = ⌊\frac{10^{2n}}{A}⌋\),然后计算\(⌊\frac{BA′}{10^{2n}}⌋\) 后调整,即得\(⌊\frac{B}{A}⌋\)

不过这里面有个\(A′\)的舍入误差问题,为了保证最后调整次数是\(O(1)\)的,那么\(A\)相对\(B\)位数不能太少,假设\(B\)\(m\)位,我们需要使得\(m ≤ 2n\),这样可以证明最后调整的时候,误差不超过\(10\)。这很简单,假如\(m > 2n\),两者同时左移若干位即可。

下面的问题就是求解\(⌊\frac{10^{2n}}{A}⌋\),设前一次迭代求解的是\(A_k^{′} = ⌊\frac{10^{2k}}{A_k}⌋\)(其中\(A_k\)\(A\)的前\(k\)位组成的数),那么这一次的迭代式为

注意这里的前\(k\)位指的是最高\(k\)位。

\[\frac{A'}{10^{2n}}=2\frac{A_k^{′}}{10^{2k}10^{n-k}}-A(\frac{A_k^{′}}{10^{2k}10^{n-k}})^2 \]

这个式子是根据\(\frac{1}{A}=\frac{A'}{10^{2n}}=\frac{A_k^{′}}{10^{2k}}\)\(\frac{1}{A}\)的牛顿迭代式\(x=2x_0-Ax_0^2\)得到的。

\[A'=2A_k^{'}10^{n-k}-\frac{A(A_k^{'})^2}{10^{2k}} \]

\[A'≈2A_k^{'}10^{n-k}-\lfloor\frac{A(A_k^{'})^2}{10^{2k}}\rfloor \]

求得\(A′\)当然还没完,这只是迭代的结果,并不是\(⌊\frac{10^{2n}}{A}⌋\)的准确值。误差分析表明,当\(k > \frac{n}{2}\)时,误差不超过\(100\)(可取\(k = ⌊\frac{n + 2}{2}⌋\)),于是还要在求得余数\(10^{2n} − AA′\)后对\(A′\)进行次数为\(O(1)\)的调整(为了降低常数,可以二分误差)。这样就迭代地求解出了\(⌊\frac{10^{2n}}{A}⌋\)

虽然跟多项式的牛顿迭代不一样这里取的是最高\(k\)位计算,但是本质是一样的,都是取估计值然后调整误差。

边界条件是\(n ≤ 2\),此时\(A\)在小整数范围内,可以直接求解结果。

注意边界条件是\(n\leq 2\)。如果只把\(n=1\)当做边界的话精度会很差。

求解出\(⌊\frac{10^{2n}}{A}⌋\)后与\(B\)相乘,再在计算余数后进行和上面类似的调整,即可求解出\(⌊\frac{B}{A}⌋\)

类似其他多项式操作,时间复杂度\(O(n\log n)\)

卡常数:乘法小范围暴力,迭代过程中没必要调整,压2位以最大限度的利用int。能做\(10^{100000}\)

#include<bits/stdc++.h>
using namespace std;

#define IN inline
#define CO const
typedef long long int64;
typedef vector<int> poly;

template<class T> IN T read(){
	T x=0,w=1;char c=getchar();
	for(;!isdigit(c);c=getchar())if(c=='-') w=-w;
	for(;isdigit(c);c=getchar()) x=x*10+c-'0';
	return x*w;
}
template<class T> IN T&read(T&x){
	return x=read<T>();
}

CO int mod=1004535809;
IN int add(int a,int b){
	return (a+=b)>=mod?a-mod:a;
}
IN int mul(int a,int b){
	return (int64)a*b%mod;
}
IN int fpow(int a,int b){
	int ans=1;
	for(;b;b>>=1,a=mul(a,a))
		if(b&1) ans=mul(ans,a);
	return ans;
}

CO int N=1<<19;
int omg[2][N],rev[N];

void NTT(poly&a,int dir){
	int lim=a.size(),len=log2(lim);
	for(int i=0;i<lim;++i){
		int r=rev[i]>>(19-len);
		if(i<r) swap(a[i],a[r]);
	}
	for(int i=1;i<lim;i<<=1)
		for(int j=0;j<lim;j+=i<<1)for(int k=0;k<i;++k){
			int t=mul(omg[dir][N/(i<<1)*k],a[j+i+k]);
			a[j+i+k]=add(a[j+k],mod-t),a[j+k]=add(a[j+k],t);
		}
	if(dir){
		int ilim=fpow(lim,mod-2);
		for(int i=0;i<lim;++i) a[i]=mul(a[i],ilim);
	}
}

CO int bas=100;

template<> poly read<poly>(){
	static char s[N];
	scanf("%s",s);
	int len=strlen(s);
	poly ans;
	for(int i=len-1;i>=0;i-=2){
		int sum=0;
		for(int j=max(i-1,0);j<=i;++j) sum=sum*10+s[j]-'0';
		ans.push_back(sum);
	}
	return ans;
}
void write(CO poly&a){
	int n=a.size()-1;
	printf("%d",a.back());
	for(int i=n-1;i>=0;--i) printf("%02d",a[i]);
}

bool operator<(CO poly&a,CO poly&b){
	if(a.size()!=b.size()) return a.size()<b.size();
	int n=a.size()-1;
	for(int i=n;i>=0;--i)if(a[i]!=b[i]) return a[i]<b[i];
	return 0;
}
bool operator<=(CO poly&a,CO poly&b){
	if(a.size()!=b.size()) return a.size()<b.size();
	int n=a.size()-1;
	for(int i=n;i>=0;--i)if(a[i]!=b[i]) return a[i]<b[i];
	return 1;
}
bool operator>(CO poly&a,CO poly&b){
	if(a.size()!=b.size()) return a.size()>b.size();
	int n=a.size()-1;
	for(int i=n;i>=0;--i)if(a[i]!=b[i]) return a[i]>b[i];
	return 0;
}
bool operator>=(CO poly&a,CO poly&b){
	return a>b or a==b;
}

poly operator<<(poly a,int m){
	int n=a.size()-1;
	a.resize(n+m+1);
	for(int i=n+m;i>=m;--i) a[i]=a[i-m];
	fill(a.begin(),a.begin()+m,0);
	return a;
}
poly operator>>(poly a,int m){
	int n=a.size()-1;
	for(int i=0;i<=n-m;++i) a[i]=a[i+m];
	a.resize(n-m+1);
	return a;
}

poly operator+(poly a,CO poly&b){
	int n=a.size()-1,m=b.size()-1;
	if(n<m) n=m,a.resize(n+1);
	for(int i=0;i<=m;++i) a[i]+=b[i];
	for(int i=0;i<n;++i)if(a[i]>=bas) ++a[i+1],a[i]-=bas;
	if(a[n]>=bas) a.push_back(1),a[n]-=bas;
	return a;
}
poly operator-(poly a,CO poly&b){ // a>=b
	int n=a.size()-1,m=b.size()-1;
	for(int i=0;i<=m;++i) a[i]-=b[i];
	for(int i=0;i<=n;++i)if(a[i]<0) --a[i+1],a[i]+=bas;
	while(a.size()>1 and !a.back()) a.pop_back();
	return a;
}
poly operator*(poly a,poly b){
	if(a==poly{0} or b==poly{0}) return poly{0};
	if(a.size()<=20 or b.size()<=20){
		int n=a.size()-1,m=b.size()-1;
		a.resize(n+m+1);
		for(int i=n;i>=0;--i){
			for(int j=m;j>=1;--j) a[i+j]+=a[i]*b[j];
			a[i]*=b[0];
		}
		for(int i=0;i<n+m;++i)if(a[i]>=bas) a[i+1]+=a[i]/bas,a[i]%=bas;
		if(a[n+m]>=bas) a.push_back(a[n+m]/bas),a[n+m]%=bas;
		return a;
	}
	int n=a.size()-1+b.size()-1,lim=1<<(int)ceil(log2(n+1));
	a.resize(lim),NTT(a,0);
	b.resize(lim),NTT(b,0);
	for(int i=0;i<lim;++i) a[i]=mul(a[i],b[i]);
	NTT(a,1),a.resize(n+1);
	for(int i=0;i<n;++i)if(a[i]>=bas) a[i+1]+=a[i]/bas,a[i]%=bas;
	if(a[n]>=bas) a.push_back(a[n]/bas),a[n]%=bas;
	return a;
}
poly operator~(poly a){
	int n=a.size();
	if(n==1){
		poly b;
		for(int x=bas*bas/a[0];x;x/=bas) b.push_back(x%bas);
		return b;
	}
	if(n==2){
		poly b;
		for(int x=bas*bas*bas*bas/(a[0]+a[1]*bas);x;x/=bas) b.push_back(x%bas);
		return b;
//		static poly v[7];
//		v[0]=a;
//		for(int i=1;i<=6;++i) v[i]=v[i-1]+v[i-1];
//		poly b(4),sum={0};
//		for(int i=3;i>=0;--i)for(int j=6;j>=0;--j)
//			if(sum+(v[j]<<i)<=(poly{1}<<2*n))
//				sum=sum+(v[j]<<i),b[i]|=1<<j;
//		while(b.size()>1 and !b.back()) b.pop_back();
//		return b;
	}
	int k=(n+2)/2;
	poly b=~poly(a.end()-k,a.end());
	b=(poly{2}*b<<(n-k))-(a*b*b>>2*k);
	return b;
}
poly operator/(poly a,poly b){
	if(a<b) return poly{0};
	int n=a.size(),m=b.size();
	if(n>2*m){
		a=a<<(n-2*m),b=b<<(n-2*m);
		n=a.size(),m=b.size();
	}
	poly c=a*~b>>2*m;
	if((c+poly{1})*b<=a){
		for(poly t=a-b*c;t>=b;t=t-b) c=c+poly{1};
	}
	else if(b*c>a){
		for(poly t=b*c-a+b-poly{1};t>=b;t=t-b) c=c-poly{1};
	}
	return c;
}

int main(){
	omg[0][0]=1,omg[0][1]=fpow(3,(mod-1)/N);
	omg[1][0]=1,omg[1][1]=fpow(omg[0][1],mod-2);
	rev[0]=0,rev[1]=1<<18;
	for(int i=2;i<N;++i){
		omg[0][i]=mul(omg[0][i-1],omg[0][1]);
		omg[1][i]=mul(omg[1][i-1],omg[1][1]);
		rev[i]=rev[i>>1]>>1|(i&1)<<18;
	}
	poly a=read<poly>(),b=read<poly>();
	write(a/b),puts("");
	return 0;
}

posted on 2020-07-17 14:58  autoint  阅读(169)  评论(0编辑  收藏

导航