P5435 基于值域预处理的快速 GCD - 数论

题意

给你 \(2n\) 个数 \(a_1,a_2,\dots,a_n;b_1,b_2,\dots,b_n\)，问 \(a,b\) 两两的 \(\gcd\)。

\(n\le 5000\)，值域 \(\le 10^6\)，时间限制 \(1\mathrm{s}\)。

题解

用 \(V\) 代表值域。

显然 \(n^2\log V\) 的做法过不了。

考虑这样一种神奇的做法：把 \(1\sim V\) 中每个整数 \(x\) 分解成三个正整数之积，且这三个数要不然 \(\le \sqrt{x}\)，要不然是质数。

具体的分解方法是，找到 \(x\) 的最小质因子 \(p\) 并获得 \(\frac{x}{p}\) 的分解，并将 \(\frac{x}{p}\) 的分解中最小的数乘上 \(p\)。

证明不会。

预处理 \(\sqrt{V}\times\sqrt{V}\) 的 \(\gcd\) 表，当分解出来的数 \(<\sqrt{V}\) 的时候直接查表，否则它一定是个质数，讨论一下就行了。

这样就可以做到 \(O(1)\) \(\gcd\) 了。

代码

#include <cstdio>
#include <cstring>
#include <cctype>
#include <cmath>
#include <algorithm>
using namespace std;
#define For(Ti,Ta,Tb) for(int Ti=(Ta);Ti<=(Tb);++Ti)
#define Dec(Ti,Ta,Tb) for(int Ti=(Ta);Ti>=(Tb);--Ti)
template<typename T> void Read(T &x){
	x=0;int _f=1;
	char ch=getchar();
	while(!isdigit(ch)) _f=(ch=='-'?-1:_f),ch=getchar();
	while(isdigit(ch)) x=x*10+(ch^48),ch=getchar();
	x=x*_f;
}
template<typename T,typename... Args> void Read(T &x,Args& ...others){
	Read(x);Read(others...);
}
typedef long long ll;
const int Inf=0x3f3f3f3f,N=5005,V=1e6+5,SV=sqrt(V),Mod=998244353;
int n,a[N],b[N];
int vis[V],prime[V],pCnt=0,k[V][3];
void Sieve(int mx){
	vis[1]=1,k[1][0]=k[1][1]=k[1][2]=1;
	For(i,2,mx){
		if(!vis[i]) prime[++pCnt]=i,k[i][0]=k[i][1]=1,k[i][2]=i;
		for(int j=1,*cur=k[i*prime[j]];1LL*i*prime[j]<=mx;++j,cur=k[i*prime[j]]){
			vis[i*prime[j]]=1;
			cur[0]=k[i][0]*prime[j],cur[1]=k[i][1],cur[2]=k[i][2];
			if(cur[0]>cur[1]) swap(cur[0],cur[1]);
			if(cur[1]>cur[2]) swap(cur[1],cur[2]);
			if(i%prime[j]==0) break;
		}
	}
}
int gcd[SV][SV];
void Pre(int mx){
	Sieve(mx);
	For(i,0,SV-1) gcd[0][i]=gcd[i][0]=i;
	For(i,1,SV-1){
		For(j,1,i){
			gcd[i][j]=gcd[j][i]=gcd[j][i%j];
		}
	}
}
int Gcd(int x,int y){
	int res=1;
	For(i,0,2){
		if(k[x][i]<SV){
			int temp=gcd[k[x][i]][y%k[x][i]];
			res*=temp,y/=temp;
		}else if(y%k[x][i]==0) res*=k[x][i],y/=k[x][i];
	}return res;
}
int main(){
	Read(n);
	int mx=0;
	For(i,1,n){
		Read(a[i]);mx=max(mx,a[i]);
	}
	For(i,1,n){
		Read(b[i]);mx=max(mx,b[i]);
	}
	Pre(mx);
	For(i,1,n){
		ll cur=i,A=0;
		For(j,1,n){
			A=(A+cur*Gcd(a[i],b[j]))%Mod;
			cur=cur*i%Mod;
		}printf("%lld\n",A);
	}
	return 0;
}