【洛谷5435】基于值域预处理的快速 GCD

题目链接

• 给定两个长度为 \(n\) 的正整数数组 \(a_{1\sim n},b_{1\sim n}\)，要求快速求出所有 \(\gcd(a_i,b_j)\)。
• \(1\le n\le 5\times10^3\)，\(1\le a_i,b_i\le10^6\)

因数分解预处理

除非一个数含有大于 \(\sqrt n\) 的质因子，否则它必然能被分解为 \(3\) 个小于等于 \(\sqrt n\) 的因数。

设 \(n\) 的最小质因子为 \(p\)（可以线性筛求出），则 \(n\) 的分解方案就是在 \(\frac np\) 的分解方案的基础上给最小的那个因数 \(\alpha\) 乘上 \(p\)。

证明：如果 \(p>\sqrt[4]n\)，显然分解出的是它的三个质因子，符合条件。否则，因为 \(\alpha\le \sqrt[3]{\frac{n}{p}}\)，所以 \(\alpha\times p\le\sqrt[3]{np^2}\le\sqrt[3]{n\cdot\sqrt n}=\sqrt n\)。

\(O(V)\) 预处理+\(O(1)\) 求 \(\gcd\)

首先预处理出 \(\sqrt V\) 范围内两两 \(\gcd\) 值 \(g_{i,j}\)。设 \(j\) 的最小质因子为 \(p\)，则 \(g_{i,j}\) 可以从 \(g_{i,j\div p}\) 转移，只需再检验 \(p\) 是否为 \(\frac{i}{g_{i,j\div p}}\) 的因子即可。

然后询问 \(\gcd\) 时就是依次求出 \(y\) 与 \(x\) 分解得到的三个因数 \(d_1,d_2,d_3\) 的 \(\gcd\)（每求出一个都要从 \(y\) 中除掉）。

如果 \(d_i\) 是大于 \(\sqrt x\) 的大质因子，只需判断 \(d_i\) 是否为 \(y\) 的因数。否则 \(\gcd(d_i,y)=\gcd(d_i,y\%d_i)\)，直接调用预处理出的值即可。

代码：\(O(V+n^2)\)

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Rg register
#define RI Rg int
#define Cn const
#define CI Cn int&
#define I inline
#define W while
#define N 5000
#define X 998244353
using namespace std;
int n,a[N+5],b[N+5];
namespace GCD
{
	#define V 1000000
	#define SV 1000
	int Pt,P[V+5],Mn[V+5],A[V+5][3],g[SV+5][SV+5];I void Init()//预处理
	{
		RI i,j;for(Mn[1]=1,i=2;i<=V;++i) for(!P[i]&&(Mn[P[++Pt]=i]=i),j=1;i*P[j]<=V;++j) if(P[i*P[j]]=1,Mn[i*P[j]]=P[j],!(i%P[j])) break;//线性筛求最小质因子
		for(A[1][0]=A[1][1]=A[1][2]=1,i=2;i<=V;++i) A[i][0]=A[j=i/Mn[i]][0]*Mn[i],A[i][1]=A[j][1],A[i][2]=A[j][2],sort(A[i],A[i]+3);//因数分解
		for(i=1;i<=SV;++i) for(j=1;j<=SV;++j) g[i][j]=i^1&&j^1?g[i][j/Mn[j]]:1,!((i/g[i][j])%Mn[j])&&(g[i][j]*=Mn[j]);//预处理SV范围内两两gcd
	}
	I int gcd(CI x,RI y)//O(1) gcd
	{
		RI o,t=1;for(RI i=0;i<=2;++i) o=y%A[x][i]?(A[x][i]<=SV?g[A[x][i]][y%A[x][i]]:1):A[x][i],y/=o,t*=o;return t;//依次求出y与x三个因数的gcd
	}
}
int main()
{
	RI i,j;for(scanf("%d",&n),i=1;i<=n;++i) scanf("%d",a+i);for(i=1;i<=n;++i) scanf("%d",b+i);
	RI t,p;for(GCD::Init(),i=1;i<=n;++i) {for(t=0,p=j=1;j<=n;++j) p=1LL*p*i%X,t=(t+1LL*p*GCD::gcd(a[i],b[j]))%X;printf("%d\n",t);}return 0;
}