bzoj 3561 DZY Loves Math VI

LINK：DZY Loves Math VI

一次询问,\(n,m\leq 500000\) 数据范围挺小。

考虑推式子 简单一推 可以推出：原式=\(\sum_{d=1}^{n}d^d\sum_{x=1}{\frac{n}{d}}x^{2d}\mu(x)\sum_{i=1}^{\frac{n}{xd}}\sum_{j=1}^{\frac{m}{xd}}i^dj^d\)

当然还可有进一步的化简为：\(\sum_{T=1}^{n}\sum_{x|T}x^{2\frac{T}{x}}\mu(x)\frac{T}{x}^{\frac{T}{x}}\sum_{i=1}^{\frac{n}{T}}i^{\frac{T}{x}}\sum_{j=1}^{\frac{n}{T}}j^{\frac{T}{x}}\)

我们发现这个是不能整除分块的 因为后面那堆东西对于d不同会变化 所以我们可以暴力枚举T 或者d 使用第一个式子或者第二个式子。

显然对于后面的那个我们枚举了d之后 变成了一个调和级数 我们前缀和一下 中间也是一个调和级数。

而第二个式子枚举T之后还需要枚举x 这个是sqrt的 没有第一个优秀 所以选择第一个式子来进行计算。

const ll MAXN=500010;
ll n,m,maxx,T,top;
ll p[MAXN],v[MAXN],mu[MAXN],g[MAXN],sum[MAXN];
inline void prepare()
{
	mu[1]=1;
	rep(2,maxx,i)
    {
        if(!v[i])
        {
            p[++top]=v[i]=i;
			mu[i]=-1;
        }
        rep(1,top,j)
        {
            if(maxx/i<p[j])break;
            ll ww=i*p[j];
            v[ww]=p[j];
            if(v[i]==p[j])break;
            mu[ww]=-mu[i];
        }
    }
}
inline ll ksm(ll b,ll p)
{
	ll cnt=1;
	while(p)
	{
		if(p&1)cnt=cnt*b%mod;
		b=b*b%mod;p=p>>1;
	}
	return cnt;
}
signed main()
{
    freopen("1.in","r",stdin);
	get(n);get(m);
	if(n>m)swap(n,m);
	ll ans=0;maxx=n;
	prepare();
	rep(1,m,i)g[i]=1;
	rep(1,n,i)
	{
		rep(1,m/i,j)g[j]=g[j]*j%mod,sum[j]=(sum[j-1]+g[j])%mod;
		ll ww=ksm(i,i);
		rep(1,n/i,j)
			ans=(ans+ww*g[j]%mod*g[j]%mod*mu[j]*sum[n/(i*j)]%mod*sum[m/(i*j)])%mod;
	}
	putl(ans);
    return 0;
}