HDU 6428 Problem C. Calculate（数论+莫比乌斯反演）

Description

求$\sum _{i=1}^A\sum _{j=1}^B\sum _{k=1}^C\phi (gcd(i,j^2,k^3))(mod{2}^{30})$

Input

第一行一整数$T$表示用例组数，每组用例输入三个整数$A,B,C(1\le T\le 10,1\le A,B,C\le 10^7)$

Output

输出答案

Sample Input

4
96 93 95
970 906 893
92460 95043 54245
9760979 8053227 7156842

Sample Output

1114536
28070648
388873924
623507672

Solution

首先有
$$\begin{array}{l}ans=\sum _{i=1}^A\sum _{d|i}\phi (d)\cdot s(i,B,C,d)\end{array}$$
其中$s(i,B,C,d)$表示满足$gcd(i,j^2,k^3)=d,1\le j\le B,1\le k\le C$的二元组$(j,k)$的个数

令$S(i,B,C,d)$表示满足$d|gcd(i,j^2,k^3),1\le j\le B,1\le k\le C$的二元组$(j,k)$的个数，那么有
$$S(i,B,C,d)=\sum _{d|n|i}s(i,B,C,k)$$
考虑求$S(i,B,C,d)$，也即分别求出求$d|j^2,1\le j\le B$的个数和$d|k^3,1\le k\le C$的个数

对$d$质因子分解有$d=p_1^{a_1}...p_m^{a_m}$，记$f(d)=p_1^{\lfloor \frac{a_1+1}{2}\rfloor }...p_m^{\lfloor \frac{a_m+1}{2}\rfloor },g(d)=p_1^{\lfloor \frac{a_1+2}{3}\rfloor }...p_m^{\lfloor \frac{a_m+2}{3}\rfloor }$，显然$f,g$均为积性函数，且$d|j^2$等价于$f(d)|j$，$d|k^3$等价于$g(d)|k$，故线性筛求出$f,g$后即可$O(1)$得到$S(i,B,C,d)=\lfloor \frac{B}{f(d)}\rfloor \cdot \lfloor \frac{C}{g(d}\rfloor$

之后由莫比乌斯反演有
$$s(i,B,C,d)=\sum _{d|n|i}\mu (\frac{n}{d})\cdot S(i,B,C,n)$$
故有
$$\begin{array}{lcl}ans& =& \sum _{i=1}^A\sum _{d|i}\phi (d)\sum _{d|n}\mu (\frac{n}{d})\cdot S(i,B,C,n)\\ & =& \sum _{i=1}^A\sum _{d|i}\phi (d)\sum _{d|n|i}\mu (\frac{n}{d})\cdot \lfloor \frac{B}{f(n)}\rfloor \cdot \lfloor \frac{C}{g(n}\rfloor \\ & =& \sum _{n=1}^A\lfloor \frac{A}{n}\rfloor \cdot \lfloor \frac{B}{f(n)}\rfloor \cdot \lfloor \frac{C}{g(n}\rfloor \sum _{d|n}\phi (d)\cdot \mu (\frac{n}{d})\end{array}$$
令$h(n)=\sum _{d|n}\phi (d)\cdot \mu (\frac{n}{d})$，由于$\phi ,\mu$都是积性函数，那么它们的狄利克雷卷积也是积性函数，故同样可以线性筛求出$h$，注意到$h(p)=\phi (1)\mu (p)+\phi (p)\mu (1)=p-2$，而由$\phi (p^c)=p^c-p^{c-1}$以及$\mu (p^c)=0,c\ge 2$可得
$$h(p^c)=\sum _{i=0}^c\phi (p^i)\mu (p^{c-i})=\phi (p^c)\mu (1)+\phi (p^{c-1})\mu (p)=p^c-2p^{c-1}+p^{c-2}$$
用线性筛$O(n)$预处理出$\phi ,\mu ,f,g,h$即可$O(A)$查询

Code

#include<cstdio>
using namespace std;
typedef long long ll;
#define maxn 10000005
#define mod (1<<30)
int mul(int x,int y)
{
	ll z=1ll*x*y;
	return z-z/mod*mod;
}
int add(int x,int y)
{
	x+=y;
	if(x>=mod)x-=mod;
	return x;
}
int prime[maxn],res,euler[maxn],f[maxn],g[maxn],h[maxn],mu[maxn],mark[maxn];
void get_euler(int n=1e7)
{
    euler[1]=mu[1]=f[1]=g[1]=h[1]=mark[1]=1;
    res=0;
    for(int i=2;i<=n;i++)
    {
        if(!euler[i])
		{
			euler[i]=i-1,prime[res++]=i,mu[i]=-1,f[i]=i,g[i]=i,h[i]=i-2,mark[i]=i;
		}
        for(int j=0;j<res&&prime[j]*i<=n;j++)
        {
            if(i%prime[j]) 
			{
				euler[prime[j]*i]=euler[i]*(prime[j]-1);
				mu[i*prime[j]]=-mu[i];
				f[i*prime[j]]=f[i]*prime[j];
				g[i*prime[j]]=g[i]*prime[j];
				mark[i*prime[j]]=prime[j];
				h[i*prime[j]]=h[i]*h[prime[j]];
			}
            else
            {
                euler[prime[j]*i]=euler[i]*prime[j];
                mu[i*prime[j]]=0;
                f[i*prime[j]]=f[i/prime[j]]*prime[j];
                if((i/prime[j])%prime[j]==0)g[i*prime[j]]=g[i/prime[j]/prime[j]]*prime[j];
                else g[i*prime[j]]=g[i];
                mark[i*prime[j]]=mark[i]*prime[j];
                if(i*prime[j]==mark[i*prime[j]])
                	h[i*prime[j]]=add(add(i*prime[j],mod-mul(2,i)),i/prime[j]);
                else
                	h[i*prime[j]]=h[i/mark[i]]*h[mark[i*prime[j]]];
                break;
            }
        }
    }
}
int main()
{
	get_euler();
	int T,A,B,C;
	scanf("%d",&T);
	while(T--)
	{
		scanf("%d%d%d",&A,&B,&C);
		int ans=0;
		for(int d=1;d<=A;d++)ans=add(ans,mul(mul(h[d],A/d),mul(B/f[d],C/g[d])));
		printf("%d\n",ans);
	}
	return 0;
}