BZOJ3561 DZY Loves Math VI

Problem

BZOJ

Solution

\[\sum_{d=1}^n \sum_{i=1}^n\sum_{j=1}^m [\gcd(i,j)=d](\frac {ij} d)^d \]

\[\sum_{d=1}^n d^d\sum_{i=1}^{n/d}\sum_{j=1}^{m/d} [\gcd(i,j)=1](ij)^d \]

考虑后面的式子

\[F(k)=\sum_{i=1}^x\sum_{j=1}^y[\gcd(i,j)=k](ij)^d \]

\[G(k)=k^{2d}\sum_{i=1}^{x/k}\sum_{j=1}^{y/k}(ij)^d \]

\[F(1)=\sum_{i=1}^x\mu(i)G(i) \]

\[\sum_{d=1}^n d^d\sum_{k=1}^{n/d}\mu(k)k^{2d}\sum_{i=1}^{n/kd}i^d\sum_{j=1}^{m/kd}j^d \]

考虑怎么暴力。枚举d再枚举k是调和级数，那么就要快速计算后面的式子，也可以通过调和级数预处理。

时间复杂度\(O(n\log n+n\ln n)\)

Code

#include <algorithm>
#include <cstdio>
using namespace std;
typedef long long ll;
const int maxn=500010,N=500000,mod=1e9+7;
template <typename Tp> inline int getmin(Tp &x,Tp y){return y<x?x=y,1:0;}
template <typename Tp> inline int getmax(Tp &x,Tp y){return y>x?x=y,1:0;}
template <typename Tp> inline void read(Tp &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
    if(ch=='-') f=1,ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
    if(f) x=-x;
}
int n,m,lim,tot,res,ans,pri[maxn],vis[maxn],mu[maxn],md[maxn],sum[maxn];
int pls(int x,int y){return x+y>=mod?x+y-mod:x+y;}
int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
int sqr(int x){return (ll)x*x%mod;}
int power(int x,int y)
{
	int res=1;
	for(;y;y>>=1,x=(ll)x*x%mod)
	  if(y&1)
	    res=(ll)res*x%mod;
	return res;
}
void init()
{
	mu[1]=md[1]=sum[1]=1;
	for(int i=2;i<=N;i++)
	{
		md[i]=1;
		if(!vis[i]) pri[++tot]=i,mu[i]=mod-1;
		for(int j=1;j<=tot&&i*pri[j]<=N;j++)
		{
			vis[i*pri[j]]=1;
			if(i%pri[j]==0){mu[i*pri[j]]=0;break;}
			mu[i*pri[j]]=dec(0,mu[i]);
		}
	}
}
int main()
{
	init();
	read(n);read(m);
	if(n>m) swap(n,m);
	for(int i=1;i<=n;i++)
	{
		lim=n/i;res=0;
		for(int j=2;i*j<=m;j++)
		{
			md[j]=(ll)md[j]*j%mod;
			sum[j]=pls(sum[j-1],md[j]);
		}
		for(int j=1;j<=lim;j++)
		  res=pls(res,(ll)mu[j]*sqr(md[j])%mod*sum[n/i/j]%mod*sum[m/i/j]%mod);
		res=(ll)res*power(i,i)%mod;
		ans=pls(ans,res);
	}
	printf("%d\n",ans);
	return 0;
}