HDU 6390 GuGuFishtion（数论+莫比乌斯反演）

Description

求$\sum\limits_{a=1}^m\sum\limits_{b=1}^n\frac{\varphi(ab)}{\varphi(a)\varphi(b)}$

Input

第一行一整数$T$表示用例组数，每组用例输入三个整数$n,m,p$，其中$p$是素数

$(1\le T\le 3,1\le n,m\le 10^6,max(m,n)<p\le 10^9+7)$

Output

输出答案，结果模$p​$

Sample Input

1
5 7 23

Sample Output

2

Solution

假设$a=p_1^{a_1}...p_r^{a_r}q_1^{b_1}...q_{s}^{b_s},b=p_1^{c_1}...p_r^{c_r}q_{s+1}^{d_1}...q_{s+t}^{d_t}$，其中$p_1,...,p_r,q_1,...,q_{s+t}$为互不相同的素数，那么有$\varphi(ab)=ab\prod\limits_{i=1}^r(1-\frac{1}{p_i})\prod\limits_{j=1}^{s+t}(1-\frac{1}{q_j})$，而$\varphi(a)\varphi(b)=ab(\prod\limits_{i=1}^r(1-\frac{1}{p_i}))^2\prod\limits_{j=1}^{s+t}(1-\frac{1}{q_j})$，故有

$$\frac{\varphi(ab)}{\varphi(a)\varphi(b)}=\frac{1}{\prod\limits_{i=1}^r(1-\frac{1}{p_i})}=\frac{gcd(a,b)}{\varphi(gcd(a,b))}$$
那么答案即为
$$\sum\limits_{i=1}^m\sum\limits_{j=1}^n\frac{gcd(i,j)}{\varphi(gcd(i,j))}=\sum\limits_{d=1}^{min(m,n)}\frac{d}{\varphi(d)}\cdot f(\lfloor\frac{m}{d}\rfloor,\lfloor\frac{n}{d}\rfloor)$$
其中$f(m,n)$表示满足$1\le a\le m,1\le b\le n,gcd(a,b)=1$的二元组$(a,b)$的对数

由莫比乌斯反演及分块加速即可$O(\sqrt{min(m,n)})$得到$f(m,n)$，注意到$1\le \varphi(d)<min(m,n)<p$，故可以线性预处理模$p$逆元，时间复杂度$O(n\sqrt{n})$

Code

#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long ll;
#define maxn 1000005
int euler[maxn],prime[maxn],res,mu[maxn],sum[maxn];
void get_euler(int n=1e6)
{
    mu[1]=sum[1]=1;
    euler[1]=1;
    res=0;
    for(int i=2;i<=n;i++)
    {
        if(!euler[i])euler[i]=i-1,prime[res++]=i,mu[i]=-1;
        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[prime[j]*i]=-mu[i];
            }
            else
            {
                euler[prime[j]*i]=euler[i]*prime[j];
                mu[prime[j]*i]=0;
                break;
            }
        }
        sum[i]=sum[i-1]+mu[i];
    }
}
int T,n,m,p;
int mul(int x,int y)
{
    ll z=1ll*x*y;
    return z-z/p*p;
}
int add(int x,int y)
{
    x+=y;
    if(x>=p)x-=p;
    return x;
}
int inv[maxn];
void init(int n)
{
    inv[0]=0;
    inv[1]=1;
    for(int i=2;i<=n;i++)inv[i]=mul(p-p/i,inv[p%i]);
}
int f(int a,int b)
{
    if(a>b)swap(a,b);
    int ans=0;
    for(int i=1,next=0;i<=a;i=next+1)
    {
        next=min(a/(a/i),b/(b/i));
        int temp=((sum[next]-sum[i-1])%p+p)%p;
        ans=add(ans,mul(mul(a/i,b/i),temp));
    }
    return ans;
}
int main()
{
    get_euler();
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d%d",&n,&m,&p);
        init(min(n,m));
        int ans=0;
        for(int d=1;d<=min(n,m);d++)
        {
            int num=f(n/d,m/d);
            ans=add(ans,mul(mul(d,num),inv[euler[d]]));
        }
        printf("%d\n",ans);
    }
    return 0;
}