Luogu P5435 基于值域预处理的快速 GCD

最近做这道题的时候被卡常了，然后突然想起来曾经偶然在陈指导的博客看到过这个\(O(1)\)做\(\gcd\)的方法

其实理解了之后还是比较简单的，以下设数的值域为\(S\)

首先我们定义对于一个数\(x\)的一种分解方式为三元组\((a,b,c)\)，满足\(a\le b\le c\and a\times b\times c=x\)，并且若\(c\)为合数则\(c\le \sqrt x\)

证明的话也很简单，考虑若\(c>\sqrt x\)，则\(a\times b<\sqrt x\)，考虑找到\(c\)的一种合法分解\(c=d\times e\)满足\(d\le e\le \sqrt x\)那么必然可以把三元组调整为\((d,a\times b,e)\)

有了这个性质我们就很好办了，考虑询问\(\gcd(x,y)\)的话我们可以分别考虑\(x\)的分解\((a,b,c)\)和\(y\)的贡献

\(a\)和\(y\)的贡献就是将答案乘上\(\gcd(a,y)\)，然后将\(y\)除以\(\gcd(a,y)\)来消去这个因子的贡献，再对\(b,c\)做同样的事即可

由于一般情况下\(a,b,c\le \sqrt x\)，因此我们可以记录下\(\sqrt S\)范围内的\(\gcd\)，这样就可以\(O(1)\)处理了，当然质数的情况可以简单讨论掉

现在就要考虑怎么\(O(S)\)求所有数的分解了，其实很简单，只要在欧拉筛的时候通过\(x\)的最小质因子\(p\)

设\(\frac{x}{p}\)的分解为\((a_0,b_0,c_0)\)，那么\(x\)的合法分解为\(\{a_0\times p,b_0,c_0\}\)的不降排序，考虑证明：

1. 当\(x\)为质数是分解显然成立
2. 当\(p\le \sqrt[4] x\)时，将\(a_0\le \sqrt[3]{\frac{x}{p}}\)代入有\(a_0\times p\le \sqrt x\)
3. 当\(p>\sqrt[4]x\)时，分情况考虑：
   • 若\(a_0=1\)，显然\(a_0\times p\le \sqrt x\)
   • 若\(a_0\ne 1\)，由于\(x\)不是质数，那么\(\frac{x}{p}\)的最小质因子\(q\)就是\(x\)的第二小质因子，一定\(\ge p\)。而\(a_0,b_0,c_0\)都为\(\frac{x}{p}\)的非\(1\)因子，故有\(p\le q\le a_0\le b_0\le c_0,p\times a_0\times b_0\times c_0>(\sqrt[4]x)^4=x\)，矛盾。故不存在这种情况

因此代码就很简单了

#include<cstdio>
#include<algorithm>
#define RI register int
#define CI const int&
using namespace std;
const int N=1000005,S=1000,mod=998244353;
int n,a[N],b[N],prm[N],cnt,G[S+5][S+5],frac[N][3]; bool vis[N];
inline void init(CI n)
{
	RI i,j; for (i=0;i<3;++i) frac[1][i]=1;
	for (i=2;i<=n;++i)
	{
		if (!vis[i]) for (prm[++cnt]=frac[i][2]=i,j=0;j<2;++j) frac[i][j]=1;
		for (j=1;j<=cnt&&i*prm[j]<=n;++j)
		{
			vis[i*prm[j]]=1; for (RI k=0;k<3;++k) frac[i*prm[j]][k]=frac[i][k];
			frac[i*prm[j]][0]*=prm[j]; sort(frac[i*prm[j]],frac[i*prm[j]]+3);
			if (i%prm[j]==0) break;
		}
	}
	for (i=1;i<=S;++i) G[i][0]=G[0][i]=i;
	for (i=1;i<=S;++i) for (j=1;j<=i;++j) G[i][j]=G[j][i]=G[i%j][j];
}
inline int gcd(int x,int y,int ret=1)
{
	for (RI i=0;i<3;++i)
	{
		int cur; if (frac[x][i]>S)
		{
			if (y%frac[x][i]) cur=1; else cur=frac[x][i];
		} else cur=G[frac[x][i]][y%frac[x][i]];
		y/=cur; ret*=cur;
	}
	return ret;
}
int main()
{
	//freopen("CODE.in","r",stdin); freopen("CODE.out","w",stdout);
	RI i,j; for (init(1000000),scanf("%d",&n),i=1;i<=n;++i) scanf("%d",&a[i]);
	for (i=1;i<=n;++i) scanf("%d",&b[i]); for (i=1;i<=n;++i)
	{
		int ans=0,cur=i; for (j=1;j<=n;++j,cur=1LL*cur*i%mod)
		(ans+=1LL*cur*gcd(a[i],b[j])%mod)%=mod; printf("%d\n",ans);
	}
	return 0;
}