調得太難受了，但我並不想多說什麼。

壓位高精

壓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) rev[i]=rev[i>>1]>>1|(i&1)<<(len-1);
	for(int i=0;i<lim;++i)if(i<rev[i]) swap(a[i],a[rev[i]]);
	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);
	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]);
	}
	poly a=read<poly>(),b=read<poly>();
	write(a/b),puts("");
	return 0;
}

LOJ164 高精度除法

壓位高精

高精度除法

題解

LOJ164 高精度除法

使用Python實現高精度除法

高精度除法

C++題解高精度除法

acwing-794. 高精度除法

c++ 超大整數除法高精度除法

高精度除法模板

【模板】高精度加法、減法、乘法、除法

內部錯誤最為致命！！！（貼一個高精度的除法）

C/C++高精度演算法的實現

如何基於Python程式碼實現高精度免費OCR工具

高精度加法、減法、乘法

9.高精度乘法

C++高精度加減乘除和求餘數

高精度減法

【演算法】【高精度】大數相關問題總結

高精度加法

高精度

AcWing 792. 高精度減法

高精度運算模板

LOJ164 高精度除法

壓位高精

高精度除法

題解

相關推薦