2018 Multi-University Training Contest 10 Calculate

题目链接：

https://vjudge.net/problem/HDU-6428

这个题题意没啥好说的，求

(sum_{i=1}^Asum_{j=1}^{B}sum_{k=1}^{C}varphi(gcd(i,j^2,k^3)))

(1le A,B,Cle 1e7)

按照这个题目的数据规模，只有线性的算法可以解决问题。

我们设(f(d,A,B,C)=sum_{i=1}^Asum_{j=1}^{B}sum_{k=1}^{C}[gcd(i,j^2,k^3)==d])

中括号里的表达式为真时，整个表达式为1，否则为0

设(F(d,A,B,C)=sum_{i=1}^Asum_{j=1}^{B}sum_{k=1}^{C}[d|gcd(i,j^2,k^3)])

即(sum_{i=1}^Asum_{j=1}^{B}sum_{k=1}^{C}[d|i\,and\,d|j^2\,and\,d|k^3])

即(sum_{i=1}^Asum_{j=1}^{B}sum_{k=1}^{C}[d|i][d|j^2][d|k^3])

即((sum_{i=1}^A[d|i])(sum_{j=1}^{B}[d|j^2])(sum_{k=1}^{C}[d|k^3]))

假设满足(d|j^2)最小的(j)设为(SQ(d))，则当仅(SQ(d)|j)时满足(d|j^2)(SQ代表平方)

假设满足(d|k^2)最小的(k)设为(CU(d))，则当仅(CU(d)|k)时满足(d|k^3)(CU代表立方)

即((lfloorfrac{A}{d} floor)(lfloorfrac{B}{SQ(d)} floor)(lfloorfrac{C}{CU(d)} floor))

如果已知质因数分解结果(n=prod_{i=1}^{m}p_i^{q_i})，有

(SQ(n)=prod_{i=1}^{m}p_i^{lceilfrac{q_i}{2} ceil},CU(n)=prod_{i=1}^{m}p_i^{lceilfrac{q_i}{3} ceil})

其中(SQ,CU)可以通过欧拉筛预处理，后续会讲

由此可知(F)可以在常数时间内计算。我们知道对(gcd()==d)的计数和对(d|gcd())的计数有莫比乌斯变换对的关系，知道一个可以求另一个。

为了在(f)和(F)之间建立关系，我们需要对已知的式子进行一些变换：

(F(d,A,B,C)=sum_{i=1}^Asum_{j=1}^{B}sum_{k=1}^{C}[d|gcd(i,j^2,k^3)]=sum_{n=1}^{A}sum_{i=1}^Asum_{j=1}^{B}sum_{k=1}^{C}[n==gcd(i,j^2,k^3)\,and\, d|n])

(=sum_{d|n\, and \, nle A}sum_{i=1}^Asum_{j=1}^{B}sum_{k=1}^{C}[n==gcd(i,j^2,k^3)]=sum_{d|n\, and \, nle A}f(d,A,B,C))

莫比乌斯变换有2种常见的形式

1.若(F(n)=sum_{d|n}f(d))

则(f(n)=sum_{d|n}mu(frac{n}{d})F(d))

2.若(F(d)=sum_{d|n \,and\,nle A}f(n))

则(f(d)=sum_{d|n \,and\,nle A}mu(frac{n}{d})F(n))

我们采取形式2，有

(f(d,A,B,C)=sum_{d|n \,and\,nle A}mu(frac{n}{d})F(n,A,B,C))

那么答案就有：

(sum_{i=1}^Asum_{j=1}^{B}sum_{k=1}^{C}varphi(gcd(i,j^2,k^3))=sum_{d=1}^Avarphi(d)f(d,A,B,C)=sum_{d=1}^Avarphi(d)sum_{d|n \,and\,nle A}mu(frac{n}{d})F(n,A,B,C))

(=sum_{n=1}^AF(n,A,B,C)sum_{d|n}varphi(d)mu(frac{n}{d}))

设(h(n)=sum_{d|n}varphi(d)mu(frac{n}{d}))，它被称为(varphi(n))和(mu(n))的狄利克雷乘积。

则答案为(sum_{n=1}^AF(n,A,B,C)h(n))

两个积性函数的狄利克雷乘积也是积性函数。

由于(varphi(n))和(mu(n))都是积性函数，所以(h(n))也是积性函数

证：若(gcd(n_1,n_2)=1)，则(h(n_1)h(n_2)=sum_{d|n_1}varphi(d)mu(frac{n_1}{d})sum_{d|n_2}varphi(d)mu(frac{n_2}{d})=sum_{d_1|n_1\,and\,d_2|n_2}varphi(d1)mu(frac{n_1}{d_1})varphi(d2)mu(frac{n_2}{d_2}))

因为(gcd(n_1,n_2)=1,d_1|n_1,d_2|n_2)，所以(gcd(d_1,d_2)=1)且(gcd(frac{n_1}{d_1},frac{n_2}{d_2})=1)

故上式为(sum_{d_1|n_1\,and\,d_2|n_2}varphi(d1d2)mu(frac{n_1}{d_1}frac{n_2}{d_2}))

即(sum d|(n_1n_2)varphi(d)mu(frac{n_1n_2}{d})=h(n_1n_2))

所以(h(n))也是积性函数

那么积性函数可以通过分解质因数的办法计算，用欧拉筛可以快速处理

由定义，(h(1)=1)

(h(p)=varphi(1)mu(p)+varphi(p)mu(1)=-1+p-1=p-2)

(h(p^q)=sum_{i=0}^{c}varphi(p^i)mu(p^{q-i})=varphi(p^{q-1})mu(p)+varphi(p^{q})mu(1)=(p-1)p^{q-2}(-1)+(p-1)p^{q-1})

(=p^{q-2}(p-1)^2,qge 2)

总结一下，只要是形如(sum_{i_1=1}^{A_1}sum_{i_2=1}^{A_2}...sum_{i_m=1}^{A_m}func(gcd(i_1^{q_1},i_2^{q_2},...,i_m^{q_m})))

其中，(func(n))是一个积性函数，比如(func(n)=n^k,varphi(n),mu(n)...)

(q_1,q_2,...,q_m,A_1,A_2,...,A_m)是正常整数，都可以用以上的方法计算，计算复杂度是(O(m*min{A_1,...,A_m}))

注：(func(p^q))的表达式不能太复杂，否则会影响计算

我们在利用欧拉筛处理时，不妨存下每个数最小的质因子和它的次数

假设我们当前处理的数字是(i)

如果(i)是一个质数，可以直接算出它的各种参数

(i)最小的质因子就是(i),次数为(1)

(SQ(i)=i)

(CU(i)=i)

(h(i)=i-2)

如果(i)是合数，那么它之前被筛过了，各种参数就算好了

用当前数(n)和质数(p)去筛选

我们保证(p)不超过(n)最小的质因子。

当(n\,mod\,p eq0)时

(n*p)最小的质因子就是(p)，次数为(1)，

(SQ(n*p)=SQ(n)*p)

(CU(n*p)=CU(n)*p)

(h(n*p)=h(n)*(p-2))

当(n\,mod\,p=0)时

(n*p)最小的质因子就是(p)，次数为(n)中(p)的次数(+1)

记这个(n*p)分解后(p)的次数为(q)

当仅(q\,mod\,2=1)时，(lceilfrac{q}{2} ceil=lceilfrac{q-1}{2} ceil+1)，有(SQ(n*p)=SQ(n)*p)

其余情况(lceilfrac{q}{2} ceil=lceilfrac{q-1}{2} ceil)，有(SQ(n*p)=SQ(n))

当仅(q\,mod\,3=1)时，(lceilfrac{q}{3} ceil=lceilfrac{q-1}{3} ceil+1)，有(CU(n*p)=CU(n)*p)

其余情况(lceilfrac{q}{3} ceil=lceilfrac{q-1}{3} ceil)，有(CU(n*p)=CU(n))

当仅(q=2)时，(gcd(n/p,p^2)=1)，有(h(n*p)=h(n/p)*h(p^2)=h(n/p)*(p-1)^2)

由于(qge 2)，其余情况一定满足(q>2)，有(h(n*p)=h(n)*p)

做此题上来用一个欧拉筛，把我们需要的这些函数都预处理出来，大约需要(5e7)级别的运算次数，

然后每次一询问都是个(O(A))的循环。由于询问10次，共计循环(1e8)次，所以要注意不要被卡常数。

看到这么大的数相加相乘，很多人会下意识地开一个(64)位整数计算，这样绝对会被卡掉。

(32)位整数的加法和乘法本质上都是对(2^{32})取模意义下的运算，所以可以开一个无符号(32)位整数，先

忽略取模，直接计算。算出来的东西自然是模(2^{32})意义下的解，之后再对(2^{30})取模即可。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<vector>
using namespace std;
const int MAXN=10000010;
bool np[MAXN];
int prime[MAXN],first[MAXN],cishu[MAXN];
int SQ[MAXN],CU[MAXN],h[MAXN];
void euler(int n)
{
	int cnt=0;
	SQ[1]=1;CU[1]=1;h[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(np[i]==0)
		{
			prime[++cnt]=i;
			first[i]=i;
			cishu[i]=1;
			SQ[i]=i;
			CU[i]=i;
			h[i]=i-2;
		}
		for(int j=1;prime[j]*i<=n&&j<=cnt;j++)
		{
			int now=prime[j]*i;
			first[now]=prime[j];
			np[now]=1;
			
			
			if(i%prime[j])
			{
				cishu[now]=1;
				SQ[now]=SQ[i]*prime[j];
				CU[now]=CU[i]*prime[j];
				h[now]=h[i]*h[prime[j]];
			}
			else
			{
				cishu[now]=cishu[i]+1;
				SQ[now]=SQ[i];
				if(cishu[now]%2==1)SQ[now]*=prime[j];
				CU[now]=CU[i];
				if(cishu[now]%3==1)CU[now]*=prime[j];
				
				if(cishu[i]==1)
				{
					h[now]=h[i/prime[j]]*(prime[j]-1)*(prime[j]-1);
				}
				else
				{
					h[now]=h[i]*prime[j];
				}
				break;
			}
		}
	}
}
const int mod=1<<30;
int main()
{
	euler(10000000);
	int T;
	scanf("%d",&T);
	while(T--)
	{
		int A,B,C;
		scanf("%d%d%d",&A,&B,&C);
		unsigned int ans=0;
		for(int i=1;i<=A;i++)
		{
			ans+=(A/i)*(B/SQ[i])*(C/CU[i])*h[i];
		}
		ans=ans%mod;
		printf("%u
",ans);
	}
}