于神之怒

题目大意

求$\sum_{i=1}^n\sum_{j=1}^mgcd(i,j)^k$
n,m,k<=5000000

繁衍

设$h[i]=i^k$，显然h可以线筛出来。
设n<=m如果不小于就调换。
设$f[d]=\sum_{i=1}^n\sum_{j=1}^m(gcd(i,j)==d)$
那么显然$f[d]=\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{m}{d}\rfloor}(gcd(i,j)==1)$
这意味着我们进行分块，每一块的f值都相同。
$ans=\sum_{d=1}^nf[d]*h[d]$
容易想到求f使用莫比乌斯反演，于是
$f[d]=\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\lfloor\frac{n}{d*i}\rfloor*\lfloor\frac{m}{d*i}\rfloor*\mu[i]$
那么同样进行分块，每一块的$\lfloor\frac{n}{d*i}\rfloor*\lfloor\frac{m}{d*i}\rfloor$一样。
所以线筛出h与$\mu$即可。

参考程序

#include<cstdio>
#include<algorithm>
#define fo(i,a,b) for(i=a;i<=b;i++)
using namespace std;
typedef long long ll;
const ll mo=1000000007;
ll mu[5000000+10],pri[1000000],h[5000000+10],sum[5000000+10],num[5000000+10];
bool bz[5000000+10];
ll i,j,k,l,r,t,n,m,nn,mm,p,top,tot,ans;
ll quicksortmi(ll x,ll y){
    if (!y) return 1;
    ll t=quicksortmi(x,y/2);
    t=t*t%mo;
    if (y%2) t=t*(x%mo)%mo;
    return t;
}
ll min(ll a,ll b){
    if (a<b) return a;else return b;
}
int main(){
    //freopen("D:/wa.txt","w",stdout);
    scanf("%lld%lld%lld",&n,&m,&p);
    if (n>m) swap(n,m);
    top=0;
    mu[1]=h[1]=1;
    fo(i,2,n){
        if (!bz[i]){
            pri[++top]=i;
            mu[i]=-1;
            h[i]=quicksortmi(i,p);
        }
        fo(j,1,top){
            if (i*pri[j]>n) break;
            h[i*pri[j]]=1LL*h[i]*h[pri[j]]%mo;
            bz[i*pri[j]]=1;
            if (i%pri[j]==0) break;
            mu[i*pri[j]]=-mu[i];
        }
    }
    sum[0]=num[0]=0;
    fo(i,1,n) num[i]=num[i-1]+mu[i],sum[i]=(sum[i-1]+h[i])%mo;
    i=1;
    while (i<=n){
        j=min(n/(n/i),m/(m/i));
        nn=n/i,mm=m/i;
        k=1;r=0;
        while (k<=nn){
            l=min(nn/(nn/k),mm/(mm/k));
            r=((nn/k)*(mm/k)%mo*(num[l]-num[k-1])%mo+r)%mo;
            k=l+1;
        }
        t=(sum[j]-sum[i-1])%mo;
        r=r*t%mo;
        ans=(ans+r)%mo;
        i=j+1;
    }
    printf("%lld\n",ans);
}