不恢复余数除法原理_大整数除法O(nlogn)

b3bbdc4d7b699eff1b938446651923ca.png

高精度加法、减法均可以在

内完成(N是位数),在FFT的帮助下高精度乘法可以做到
(或许是
?),那么加减乘除就剩下了除法,我们可以试着把除法的复杂度降下来。

Naive的想法1:

既然直接列竖式做除法是

的,那么我们是不是可以二分答案呢?

仔细想想,二分一个

,乘法一个
,所以复杂度是
,稳稳的。

哎呀,二分的

是与位数同级的(即
),到头来还是
的。

Naive的想法2:

既然多项式除法可以做到

,那么是不是整数除法也可以看做是多项式除法呢?

不行,两者定义有差异(多项式除法要求余数度数小于被除数而整数要求余数数值上小于被除数)。

正确做法:

二分的问题出在哪里?要求答案达到

位精度需要二分
次,在大数据下即使是
也会很大。比二分更快的是牛顿迭代法,同样使答案达到
位精度只需要迭代
次。

我们要求

,可以先求
再乘上
,这样化出来的式子不包含除法:

时间复杂度:

初始时有效位数为1,之后每迭代一次位数翻倍,乘法使用FFT优化后,时间复杂度满足

(n为有效位数)。

实现细节:

1.为了方便可以将整数转换成多项式

2.类似多项式求逆,每一次计算时只需要取

的前
位参加计算,同样答案需要舍弃后面的位数。

3.迭代初值在

收敛,那么初值可以设置为答案的第一位(例如
时初值设为
,
时,初值设为
)。
例:求

转换问题变成求
。先求
,初值设为
,

第一次迭代:

第二次迭代:

第三次迭代:

第四次迭代:

第五次迭代

需要注意的是最后的答案只有一半的位数是正确的,计算答案得到

下面是NTT实现的高精度除法(答案向下取整)。

#include<algorithm> 
#include<iostream>
#include<cstring>
#include<cstdio>
#include<vector>
#define loop(n,i) for(register int i=1;i<=(n);++i)
#define MAX 500009
#define zxc(x) cerr<<(#x)<<'='<<(x)<<'n'
#define zxcv(x) cerr<<(#x)<<'='<<(x)<<','
#define int long long
#define P 998244353
using namespace std;
inline int icin(){
	char c=getchar();int s=0;bool sign=0;
	while(!isdigit(c)&&c^'-')c=getchar();
	if(c=='-')c=getchar(),sign=1;
	while(isdigit(c))s=(s<<1)+(s<<3)+c-'0',c=getchar();
	return sign?-s:s;
}
typedef vector<int> Poly;
inline int Quick(int a,int m){int ans=1;for(;m;m>>=1,a=a*a%P) if(m&1) ans=ans*a%P;return ans;}
int R[MAX];
inline void NTT(Poly& a,int x,int type){
	int n=1<<x;
	for(register int i=0;i<n;++i) if((R[i]=R[i>>1]>>1|(i&1)<<x-1)>i) swap(a[i],a[R[i]]);
	for(register int s=1;s<=x;++s){
		int len=1<<s,mid=len>>1;
		int w=Quick(3,type*(P-1)/len+P-1);
		for(register int k=0;k<n;k+=len){
			int d=1;
			for(register int j=0;j<mid;++j){
				int u=a[j+k],v=a[j+k+mid]*d%P;
				a[j+k]=(u+v)%P,a[j+k+mid]=(u-v+P)%P;
				d=d*w%P;
			}
		}
	}
	if(type==-1) for(register int i=0,inv=Quick(n,P-2);i<n;++i) a[i]=a[i]*inv%P;
}
ostream& operator <<(ostream& out,Poly x){for(register int i=0;i<x.size();++i) out<<x[i]<<' ';return out;}
inline int ceilog(int x){int ans=0;while(1<<ans<x)ans++;return ans;}
Poly operator *(Poly A,Poly B){
	int deg=A.size()+B.size()-1;
	int x=ceilog(deg);
	A.resize(1<<x),B.resize(1<<x);
	NTT(A,x,1),NTT(B,x,1);
	for(register int i=0;i<1<<x;++i) A[i]=A[i]*B[i]%P;
	NTT(A,x,-1);
	A.resize(deg);
	for(register int i=deg-1;i;--i){
		A[i-1]+=A[i]/10;
		A[i]%=10;
	}
	return A;
}

inline Poly Inverse(Poly A){
	Poly B,C;B.resize(2);B[1]=100/(A[0]*10+(A.size()>1?A[1]:0));
	int x=ceilog(A.size())+1;
	for(register int s=1;s<=10;++s){//迭代次数 
		C.resize(1<<s),B.resize(1<<s);
		for(register int i=0;i<min(1ll<<s,(int)A.size());++i) C[i]=A[i];
		for(register int i=min(1ll<<s,(int)A.size());i<1<<s;++i) C[i]=0;
		C=B*C;
		for(register int i=1;i<C.size();++i) C[i]=-C[i];
		C[0]=2-C[0];
		for(register int i=C.size()-1;i;--i){
			C[i-1]+=(C[i]-9)/10;
			C[i]=(C[i]+10)%10;
		}
		B=B*C;
		B.resize(1<<s);
	}
	return B;
}
Poly A,B;
int rest;
main(){
	char c=getchar();while(!isdigit(c)) c=getchar();
	while(isdigit(c)) A.push_back(c-'0'),c=getchar();
	rest+=A.size();
	c=getchar();while(!isdigit(c)) c=getchar();
	while(isdigit(c)) B.push_back(c-'0'),c=getchar();
	rest-=B.size();
	A=A*Inverse(B);
	for(register int i=0;i<=rest;++i) if(!(i==0&&A[i]==0)) cout<<A[i];
	cout<<'n';
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值