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;
}
静渊以有谋，疏通而知事。
