[bzoj3561] DZY Loves Math VI

题目大意

给定n,m，求$\sum_{i=1}^n \sum_{j=1}^m [i,j]^{(i,j)}$模$10^9+7$的值。

n,m≤500000

分析

这题看起来不多人过，但是感觉挺简单的。。。
令n≤m
首先枚举gcd的值，然后得到：

$$Ans=\sum_{d=1}^n \sum_{i=1}^{\lfloor\frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor\frac{m}{d} \rfloor}(ijd)^d*[(i,j)=1]$$
然后直接上莫比乌斯反演，得到：
$$Ans=\sum_{d=1}^n d^d \sum_{d'=1}^{\lfloor\frac{n}{d} \rfloor}\mu(d') [(\sum_{i=1}^{\lfloor\frac{n}{dd'} \rfloor}(id')^d)*(\sum_{j=1}^{\lfloor\frac{m}{dd'} \rfloor}(jd')^d)]$$
还要继续推吗？
数据范围只是500000，上面的式子是可以mlogm算的。而且继续我也不会推。

#include <cstdio>
#include <cstring>
#include <algorithm>

using namespace std;

const int N=5e5,mo=1e9+7;

typedef long long LL;

int T,mu[N+5],tot,p[N],f[N+5],n,m,ans,s[N+5],t[N],s1,s2;

bool bz[N+5];

int quick(int x,int y)
{
    if (!y) return 1;
    int ans=quick(x,y>>1);
    ans=(LL)ans*ans%mo;
    if (y&1) ans=(LL)ans*x%mo;
    return ans;
}

int main()
{
    scanf("%d%d",&n,&m);
    if (n>m) n^=m^=n^=m;
    mu[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!bz[i])
        {
            mu[i]=-1;
            p[tot++]=i;
        }
        for (int j=0;j<tot && i*p[j]<=n;j++)
        {
            int I=i*p[j];
            bz[I]=1;
            if (i%p[j]==0)
            {
                mu[I]=0; break;
            }
            mu[I]=-mu[i];
        }
    }
    for (int i=1;i<=m;i++) t[i]=1;
    for (int d=1;d<=n;d++)
    {
        for (int i=1;i*d<=m;i++) t[i]=(LL)t[i]*i%mo,s[i]=(s[i-1]+t[i])%mo;
        int sum=0;
        for (int i=1;i*d<=n;i++) sum=(sum+(LL)s[n/(i*d)]*s[m/(i*d)]*mu[i]%mo*t[i]%mo*t[i])%mo;
        ans=(ans+(LL)sum*quick(d,d))%mo;
    }
    ans=(ans+mo)%mo;
    printf("%d\n",ans);
    return 0;
}