标签：LOJ164 inter return bas 高精度 int poly frac 除法
调得太难受了，但我并不想多说什么。
压位高精
压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,inter,return,bas,高精度,int,poly,frac,除法 来源： https://www.cnblogs.com/autoint/p/13330002.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。