BZOJ4454. C Language Practice GCD的O(n)预处理O(1)求法

题意：

给定一个长度为\(n,m\)的序列\(a,b\)，求\(\sum_{i=0}^{n-1}\sum_{j=0}^{m-1}gcd(a_i,b_j)xor\;i\;xor \;j\)的值，结果对\(2^{32}\)取模,多组数据

范围&性质：\(1\le n,m\le 2*10^3,1\le a_i,b_j\le 10^6,1\le t\le 85\)

分析：

一眼暴力做法：复杂度 \(\omicron(tn^2log)\approx 8*10^9\) 铁定被卡

我们分析可得\(n^2\)的枚举无法优化，\(t\)无法改变，那么唯一能降低复杂度的地方就是GCD的\(log\)，而且本题的数据范围里特意给出了\(a_i，b_j\)的大小，不超过\(10^6\)，显然我们要通过一些方法将GCD的复杂度降到\(\omicron(1)\)

---

以下方法出自Claris的论文：

1. 首先我们记\(\sqrt[2]{n}\)的大小为\(m\)，对于\(1\le i,j\le m\)的所有数求一遍GCD，复杂度为\(\omicron(m*m)=\omicron(n)\)

2. 然后我们对于任意一个不超过\(n\)的数\(x\)，一定存在一种方法，使得\(x\)分解成三个数\(a,b,c\)，满足\(a,b,c\)要么小于等于\(m\),要么是质数（证明详见原作者论文）

3. 那么\(x,y\)的GCD可以看成\(a,b,c\)与\(y\)的合并

   若\(a\)为质数且y为a的倍数,\(y/=a,ans*=a\)

   若\(a\le m\)，\(gcd(y,a)=gcd（a,y\;mod\;a)\)而\(gcd（a,y\;mod\;a)\)已经预处理过了

像这样就可以做到\(\omicron(n)\)预处理\(\omicron(1)\)查询，注意这道题目卡内存，只能精确取\(m=\sqrt[2]{n}\)，实践中可以把m取得大一点。

代码：

#include<bits/stdc++.h>

using namespace std;

namespace zzc
{
	const int maxn =1e6+5;
	int t,n,m,cnt=0;
	int g[1005][1005],p[100005],num[maxn],f[maxn][3];
	int a[2005],b[2005];
	bool vis[maxn];
	
	int getgcd(int x,int y)
	{
		return (x&&y)?g[y][x%y]:x|y;
	}
	
	void init()
	{
		num[1]=1;
		for(int i=2;i<=1000000;i++)
		{
			if(!vis[i])
			{
				p[++cnt]=i;
				num[i]=i;
			}
			for(int j=1;j<=cnt&&i*p[j]<=1000000;j++)
			{
				vis[i*p[j]]=true;
				num[i*p[j]]=p[j];
				if(i%p[j]==0) break;
			}
		}
	}
	
	int gcd(int x,int y)
	{
		if(x<=1000&&y<=1000) return g[x][y];
		else if(!x||!y) return x|y;
		int d=1;
		for(int i=0;i<3;i++)
		{
			if(f[x][i]==1) continue;
			int tmp=f[x][i];
			if(num[tmp]==tmp)
			{
				if(y%tmp==0)
				{
					y/=tmp;
					d*=tmp;
				}
			}
			else
			{
				int t=g[tmp][y%tmp]; 
				y/=t;
				d*=t;
			}
		} 
		return d;
	}
	
	void work()
	{
		for(int i=0;i<=1000;i++)
		{
			for(int j=0;j<=i;j++)
			{
				g[i][j]=g[j][i]=getgcd(i,j);
			}
		}
		init();
		f[1][0]=f[1][1]=f[1][2]=1;
		for(int i=2;i<=1000000;i++)
		{
			memcpy(f[i],f[i/num[i]],sizeof(f[i]));
			if(f[i][0]*num[i]<=1000) f[i][0]*=num[i];
			else if(f[i][1]*num[i]<=1000) f[i][1]*=num[i];
			else f[i][2]*=num[i];
		}
		scanf("%d",&t);
		while(t--)
		{
			scanf("%d%d",&n,&m);
			for(int i=0;i<n;i++)
			{
				scanf("%d",&a[i]);
			}
			for(int j=0;j<m;j++)
			{
				scanf("%d",&b[j]);
			}
			unsigned int ans=0;
			for(int i=0;i<n;i++)
			{
				for(int j=0;j<m;j++)
				{
					ans+=gcd(a[i],b[j])^i^j;
				}
			}
			printf("%u\n",ans);
		}
		
	}
	
}

int main()
{
	zzc::work();
	return 0;
}