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

题目

题目链接：https://www.luogu.com.cn/problem/P5435
给定 (n) 个正整数 (a_1,a_2,dots,a_n)，再给定 (n) 个正整数(b_1,b_2,dots,b_n)，你需要对每对 ((i,j)) 求出 (a_i) 与 (b_j) 的最大公因数。
不难发现你的输出应有 (n^2) 个正整数。为了减少输出对程序的运行效率的影响，你只需要输出 (n) 行，每行一个整数 (A_i)。
其中对于 (iin[1,n])，(A_i=sum_{j=1}^{n}i^jgcd(a_i,b_j))。由于答案可能过大，你只需要输出模 (998,244,353) 后的结果即可。
(nleq 5000)，(1leq a_ileq 10^6)。

思路

一种挺有趣的科技。
对于任意一个数 (x)，肯定有一种方法可以把 (x) 分解成 (x_1 imes x_2 imes x_3)，且 (x_1,x_2,x_3) 都 (leq sqrt{x}) 或是质数。
具体的构造方法，我们可以线性筛出值域内所有数的最小质因子。对于一个数 (p)，如果它是质数，显然可以分解为 (1 imes 1 imes p)；如果它不是质数，设它的最小质因子是 (d)，且 (q=frac{p}{d}) 分解为了 (q_1 imes q_2 imes q_3)（(q_1leq q_2leq q_3)），那么可以把 (p) 分解为 (q_1 imes d,q_2,q_3)。
证明：

• 如果 (dleq sqrt[4]{p})，因为 (q_1leq sqrt[3]{frac{p}{d}})，所以 (d imes q_1leq d imes sqrt[3]{frac{p}{d}}=p^{frac{1}{3}} imes d^{frac{2}{3}})。而 (sqrt{p}div p^{frac{1}{3}}=p^{frac{1}{6}})，(d^{frac{2}{3}}leq left(sqrt[4]{p} ight)^{frac{2}{3}}leq p^{frac{1}{6}})。故 (d imes q_1leq sqrt{p})。
• 如果 (d>sqrt[4]{p})
  • 如果 (q_1=1)，显然 (1 imes dleq sqrt{p})。
  • 如果 (q_1>1)，则 (dleq q_1leq q_2leq q_3)，则 (d imes q_1 imes q_2 imes q_3>left(sqrt[4]{p} ight)^4=p)，矛盾。

这样把值域内所有的数分解好后，预处理出一个 (sqrt{v} imes sqrt{v}) 的 (gcd) 数组（(v) 是值域）。对于每一次询问 (a,b)，把 (a) 分解成 (a_1 imes a_2 imes a_3)。

• 如果 (a_1) 是质数：如果 (b) 是 (a_1) 的倍数返回 (a_1)，否则返回 (1)。
• 如果 (a_1) 不是质数，(gcd(a_1,b)=gcd(a_1,bmod a_1))，此时两个数都不超过 (sqrt{v})，直接返回预处理的值。

对 (a_1,a_2,a_3) 分别做一次乘起来，就是 (gcd(a,b)) 了。注意每次 (b) 需要除掉返回的 (gcd)。
时间复杂度 (O(v+n^2))。

代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

const int N=5010,M=1000010,MM=1010,MOD=998244353;
int n,m,a[N],b[N],prm[M],v[M],f[M][3],g[MM][MM];

void findprm(int n)
{
	v[1]=1;
	for (int i=2;i<=n;i++)
	{
		if (!v[i]) prm[++m]=i,v[i]=i;
		for (int j=1;j<=m;j++)
		{
			if (i>n/prm[j] || prm[j]>v[i]) break;
			v[i*prm[j]]=prm[j];
		}
	}
}

void prework()
{
	for (int i=1;i<MM;i++) g[0][i]=g[i][0]=i;
	for (int i=1;i<MM;i++)
		for (int j=1;j<=i;j++)
			g[i][j]=g[j][i]=g[j][i%j];
	for (int i=1;i<M;i++)
		if (v[i]==i)
			f[i][0]=1,f[i][1]=1,f[i][2]=i;
		else
		{
			f[i][0]=f[i/v[i]][0]*v[i]; f[i][1]=f[i/v[i]][1]; f[i][2]=f[i/v[i]][2];
			if (f[i][0]>f[i][1]) swap(f[i][0],f[i][1]);
			if (f[i][1]>f[i][2]) swap(f[i][1],f[i][2]);
		}
}

int gcd(int x,int &y)
{
	int d=(x>MM)?((y%x)?1:x):g[x][y%x];
	y/=d;
	return d;
}

int main()
{
	findprm(M-1);
	prework();
	scanf("%d",&n);
	for (int i=1;i<=n;i++) scanf("%d",&a[i]);
	for (int i=1;i<=n;i++) scanf("%d",&b[i]);
	for (int i=1;i<=n;i++)
	{
		ll ans=0,res=1;
		for (int j=1,k;j<=n;j++)
		{
			res=res*i%MOD; k=b[j];
			ans=(ans+res*gcd(f[a[i]][0],k)*gcd(f[a[i]][1],k)*gcd(f[a[i]][2],k))%MOD;
		}
		cout<<ans<<"
";
	}
	return 0;
}