调得太难受了，但我并不想多说什么。
压位高精
压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 floor)。
对于100%的数据，(1leq bleq 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}}=2frac{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}-lfloorfrac{A(A_k^{'})^2}{10^{2k}}
floor
]
求得(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)在小整数范围内，可以直接求解结果。
注意边界条件是(nleq 2)。如果只把(n=1)当做边界的话精度会很差。
求解出(⌊frac{10^{2n}}{A}⌋)后与(B)相乘，再在计算余数后进行和上面类似的调整，即可求解出(⌊frac{B}{A}⌋)。
类似其他多项式操作，时间复杂度(O(nlog 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;
}