BZOJ4407：于神之怒加强版（线性筛）

题面
题意：对于给定的k，求

$$\sum_{i=1}^n\sum_{j=1}^mgcd(i,j)^k$$
n，m，k≤5e6，2000组数据。

根据路人试子的推导（倒）方法，我们枚举gcd，得
$d^k$完全积性，可以先处理质数的结果，再算出全部。
且由于这个完全积性，它卷个$\mu$也可以用同样的方法算出来，然后对于每个询问分块就可以了。

#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <cstdio>
#include <cstdlib>
#include <cstring>

using namespace std;
#define mmst(a, b) memset(a, b, sizeof(a))
#define mmcp(a, b) memcpy(a, b, sizeof(b))

typedef long long LL;

const int N=5005000,nn=5000000;
const LL mo=1e9+7;

int T,k;
int miu[N],u[N],prime[N],num;
LL n,m,h[N],dk[N];
bool b[N];

LL cheng(LL a,LL b)
{
    LL res=1;
    for(;b;b>>=1,a=a*a%mo)
    if(b&1)
    res=res*a%mo;
    return res;
}

int main()
{
    cin>>T>>k;

    h[1]=dk[1]=miu[1]=1;
    for(int i=2;i<=nn;i++)
    {
        if(!b[i])
        {
            u[i]=i;
            miu[i]=-1;
            dk[i]=cheng(i,k);
            h[i]=(dk[i]-1+mo)%mo;
            prime[++num]=i;

        }
        for(int j=1;j<=num&&prime[j]*i<=nn;j++)
        {
            int hy=i*prime[j];
            u[hy]=u[i];
            b[hy]=1;
            dk[hy]=dk[i]*dk[prime[j]]%mo;   
            if(i%prime[j]==0)
            {
                miu[hy]=0;
                h[hy]=h[i]*dk[prime[j]]%mo;
                break;
            }
            miu[hy]=-miu[i];
            h[hy]=h[i]*h[prime[j]]%mo;
        }
    }

    for(int i=2;i<=nn;i++)
    h[i]=(h[i-1]+h[i])%mo;

    while(T--)
    {
        scanf("%lld%lld",&n,&m);
        if(n>m)
        swap(n,m);
        LL ans=0;
        for(LL last,i=1;i<=n;i=last+1)
        {
            last=min(n/(n/i),m/(m/i));
            ans=(ans+(n/i)*(m/i)%mo*(h[last]-h[i-1]+mo)%mo )%mo;
        }
        printf("%lld\n",ans);
    }

    return 0;
}