BZOJ3561 DZY Loves Math VI

原题链接：https://www.lydsy.com/JudgeOnline/problem.php?id=3561

DZY Loves Math VI

Description

给定正整数n,m。求

$$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)^{gcd(i,j)}$$

Input

一行两个整数n,m。

Output

一个整数，为答案模1000000007后的值。

Sample Input

5 4

Sample Output

424

HINT

数据规模：

1<=n,m<=500000，共有3组数据。

题解

$$\begin{align*} &\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)^{gcd(i,j)}\\ &=\sum_{i=1}^{n}\sum_{j=1}^{m}(\frac{ij}{gcd(i,j)})^{gcd(i,j)}\\ &=\sum_{d=1}^{min(n,m)}\sum_{i=1}^{n}\sum_{j=1}^{m}(\frac{ij}{d})^{d}[d=gcd(i,j)]\\ &=\sum_{d=1}^{min(n,m)}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}(\frac{d^2ij}{d})^{d}[1=gcd(i,j)]\\ &=\sum_{d=1}^{min(n,m)}d^d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}(ij)^{d}\sum_{d'|gcd(i,j)}\mu(d')\\ &=\sum_{d=1}^{min(n,m)}d^d\sum_{d'=1}^{\lfloor\frac{min(n,m)}{d}\rfloor}\mu(d')\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}(ij)^{d}[d'|gcd(i,j)]\\ &=\sum_{d=1}^{min(n,m)}d^d\sum_{d'=1}^{\lfloor\frac{min(n,m)}{d}\rfloor}\mu(d')d'^{2d}\sum_{i=1}^{\lfloor\frac{n}{dd'}\rfloor}i^d\sum_{j=1}^{\lfloor\frac{m}{dd'}\rfloor}j^{d}\\ \end{align*}$$

因为$dd'$后面还跟了一个$i^d$之类的奇奇怪怪的东西，所以用$T$来替换并无卵用，但是数据范围只有$5e4$啊，如果直接枚举的话复杂度为$\frac{n}{1}+\frac{n}{2}+\frac{n}{3}+\cdots+\frac{n}{n-1}+\frac{n}{n}\approx n\ ln(n)$，稳了稳了。

代码

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int M=5e5+5,mod=1000000007;
int p[M/3],miu[M],n,m,nd,md,dd;
ll ans,tmp,sum[M],val[M];
bool isp[M];
ll power(ll x,ll p){ll r=1;for(;p;p>>=1,x=x*x%mod)if(p&1)r=r*x%mod;return r;}
void pre()
{
    miu[1]=isp[1]=1;
    for(int i=2,t;i<M;++i)
    {
        if(!isp[i])p[++p[0]]=i,miu[i]=-1;
        for(int j=1;j<=p[0],i*p[j]<M;++j){isp[i*p[j]]=1;if(i%p[j]==0)break;miu[i*p[j]]=-miu[i];}
    }
}
void in(){scanf("%d%d",&n,&m);}
void ac()
{
    if(n>m)swap(n,m);pre();
    for(int i=1;i<=m;++i)val[i]=1;
    for(int i=1,j;i<=n;++i)
    {
        nd=n/i,md=m/i,tmp=0;
        for(j=1;j<=md;++j)val[j]=val[j]*j%mod,sum[j]=(sum[j-1]+val[j])%mod;
        for(j=1;j<=nd;++j)(tmp+=miu[j]*val[j]%mod*val[j]%mod*sum[nd/j]%mod*sum[md/j]%mod)%=mod;
        (ans+=tmp*power(i,i)%mod)%=mod;
    }
    printf("%lld",(ans+mod)%mod);
}
int main(){in();ac();}