题目链接：
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);
}
}