【BZOJ3561】DZY Loves Math VI （数论）

【BZOJ3561】DZY Loves Math VI （数论）

题面

BZOJ

题解

\[\begin{aligned} ans&=\sum_{i=1}^n\sum_{j=1}^m\sum_{d=1}^n[gcd(i,j)=d](\frac{ij}{d})^d\\ &=\sum_{d=1}^nd^d\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)=1]i^dj^d\\ &=\sum_{d=1}^nd^d\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}i^dj^d\sum_{x|gcd(i,j)}\mu(x)\\ &=\sum_{d=1}^nd^d\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}i^dj^d\sum_{x|i,x|j}\mu(x)\\ &=\sum_{d=1}^nd^d\sum_{x=1}^{n/d}\mu(x)\sum_{i=1}^{n/xd}\sum_{j=1}^{m/xd}(ix)^d(jx)^d\\ &=\sum_{d=1}^nd^d\sum_{x=1}^{n/d}\mu(x)x^{2d}\sum_{i=1}^{n/xd}\sum_{j=1}^{m/xd}i^dj^d\\ \end{aligned}\]

然后发现\(\sum_i i^d\)不会算，实际上枚举\(d\)的时候就大力预处理一次，这样子的预处理的复杂度是调和级数的。
然后整个式子都调和级数的爆算就完了。。

#include<iostream>
#include<cstdio>
using namespace std;
#define MOD 1000000007
#define MAX 500500
int mu[MAX],pri[MAX],tot;
bool zs[MAX];
int n,m,ans,v[MAX],s[MAX],x[MAX],D[MAX];
int fpow(int a,int b){int s=1;while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}return s;}
void pre(int n)
{
	mu[1]=1;
	for(int i=2;i<=n;++i)
	{
		if(!zs[i])pri[++tot]=i,mu[i]=MOD-1;
		for(int j=1;j<=tot&&i*pri[j]<=n;++j)
		{
			zs[i*pri[j]]=true;
			if(i%pri[j]==0){mu[i*pri[j]]=0;break;}
			mu[i*pri[j]]=MOD-mu[i];
		}
	}
}
int main()
{
	scanf("%d%d",&n,&m);if(n>m)n^=m,m^=n,n^=m;pre(n);
	for(int i=1;i<=m;++i)v[i]=x[i]=1;
	for(int i=1;i<=n;++i)D[i]=fpow(i,i);
	for(int d=1;d<=n;++d)
	{
		for(int i=1;i<=m/d;++i)v[i]=1ll*v[i]*i%MOD;
		for(int i=1;i<=m/d;++i)s[i]=(s[i-1]+v[i])%MOD;
		for(int i=1;i<=n/d;++i)x[i]=1ll*x[i]*i%MOD*i%MOD;
		for(int i=1;i<=n/d;++i)ans=(ans+1ll*D[d]*mu[i]%MOD*x[i]%MOD*s[n/d/i]%MOD*s[m/d/i])%MOD;
	}
	printf("%d\n",ans);
	return 0;
}