[BZOJ3561] DZY Loves Math VI

Description

给定正整数n,m。求

Input

一行两个整数n,m。

Output

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

Sample Input

5 4

Sample Output

424

Solution

推下式子，莫比乌斯反演一波：

\[\begin{align} ans=&\sum_{i=1}^n\sum_{j=1}^{m}(\frac{ij}{\gcd(i,j)})^{\gcd(i,j)}\\ =&\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m}(\frac{ij}{d})^d[gcd(i,j)=d]\\ =&\sum_{d=1}^nd^d\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}(ij)^d\sum_{t|i,t|j}\mu(t)\\ =&\sum_{d=1}^nd^d\sum_{t=1}^{n/d}\mu(t)t^{2d}\sum_{i=1}^{n/dt}i^d\sum_{j=1}^{m/dt}j^d\\ \end{align} \]

注意上面推导默认\(n\leqslant m\)。

然后一路暴力算就好了，根据调和级数，复杂度为\(O(n\log n)\)。

注意算幂的时候写的好一点，不要每次快速幂，否则复杂度多个\(\log\)就\(T\)了。

#include<bits/stdc++.h>
using namespace std;
 
void read(int &x) {
    x=0;int f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}
 
void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('\n');}

#define lf double
#define ll long long 

const int maxn = 5e5+10;
const int inf = 1e9;
const lf eps = 1e-8;
const int mod = 1e9+7;

int n,m,vis[maxn],pri[maxn],tot,mu[maxn],sum[maxn],pw[maxn];

void sieve() {
	mu[1]=1;
	for(int i=2;i<maxn;i++) {
		if(!vis[i]) pri[++tot]=i,mu[i]=-1;
		for(int j=1;j<=tot&&i*pri[j]<maxn;j++) {
			vis[i*pri[j]]=1;
			if(i%pri[j]==0) break;
			mu[i*pri[j]]=-mu[i];
		}
	}
}

int qpow(int a,int x) {
	int res=1;
	for(;x;x>>=1,a=1ll*a*a%mod) if(x&1) res=1ll*res*a%mod;
	return res;
}

int main() {
	sieve();
	read(n),read(m);int ans=0;
	if(n>m) swap(n,m);
	for(int i=1;i<=m;i++) pw[i]=1;
	for(int d=1;d<=n;d++) {
		int res=0;sum[0]=0;
		for(int t=1;t<=m/d;t++) pw[t]=1ll*pw[t]*t%mod;
		for(int t=1;t<=m/d;t++) sum[t]=(sum[t-1]+pw[t])%mod;
		for(int t=1;t<=n/d;t++) res=(res+1ll*mu[t]*pw[t]*pw[t]%mod*sum[n/d/t]%mod*sum[m/d/t]%mod)%mod;
		res=1ll*res*qpow(d,d)%mod;
		ans=(ans+res)%mod;
	}write((ans+mod)%mod);
	return 0;
}