大致题意： 求(sum_{i=1}^nsum_{j=1}^mgcd(i,j)^k)。
前言
(Jan 28th)刷题计划（3/6），算法标签：莫比乌斯反演&杜教筛。
学了一个新技巧：线性筛筛自定义积性函数。
莫比乌斯反演
首先按照套路去枚举(gcd)得到：
[sum_{d=1}^{min(n,m)}d^ksum_{i=1}^{lfloorfrac nd
floor}sum_{j=1}^{lfloorfrac md
floor}[gcd(i,j)==1]
]
莫比乌斯反演，得到：
[sum_{d=1}^{min(n,m)}d^ksum_{p=1}^{min(lfloorfrac nd
floor,lfloorfrac md
floor)}mu(p)lfloorfrac n{dp}
floorlfloorfrac m{dp}
floor
]
考虑枚举(D=dp)，得到：
[sum_{D=1}^{min(n,m)}lfloorfrac nD
floorlfloorfrac mD
floorsum_{d|D}d^kmu(frac Dd)
]
设(f(D)=sum_{d|D}d^kmu(frac Dd))，则原式等同于：
[sum_{D=1}^{min(n,m)}lfloorfrac nD
floorlfloorfrac mD
floor f(D)
]
则只要求出(f(D))，就可以通过除法分块轻松解决此题。
于是关键问题来了：如何求出(f(D))。
线性筛筛积性函数
显然(f(n))是一个积性函数。
如果我们设(n=prod_{i=1}^tP_i^{a_i})（(P)为质数），那么：
[f(n)=prod_{i=1}^tf(P_i^{a_i})
]
考虑化简(f(P_i^{a_i}))：
[f(P_i^{a_i})=sum_{d|P_i^{a_i}}d^kmu(frac Dd)=1^k imesmu(P_i^{a_i})+P_i^k imesmu(P_i^{a_i-1})+...+P_i^{a_i imes k} imesmu(1)
]
根据(mu)的定义可知，(mu(P_i^{a_i}),mu(P_i^{a_i-1}),...,mu(P_i^2))都等于(0)，也就是说，我们只需要考虑这个式子的最后两项。即：
[f(P_i^{a_i})=P_i^{(a_i-1) imes k} imesmu(P_i)+P_i^{a_i imes k} imesmu(1)=-P_i^{(a_i-1) imes k}+P_i^{a_i imes k}=P_i^{(a_i-1) imes k}(P_i^k-1)
]
也就是说，有：
[f(n)=prod_{i=1}^tP_i^{(a_i-1) imes k}(P_i^k-1)
]
然后就很好做了。
以线性筛为例，我们可以得到边界状态为：
[f(n)=egin{cases}1&(n=1)\p^k-1&(n is prime)end{cases}
]
然后考虑筛的过程中的筛法为：
[f(i imes P_j)=egin{cases}f(i) imes f(P_j)&(P_j
ot{|}i)\f(i) imes P_i^k&(P_j|i)end{cases}
]
这可以根据前面推出的式子自行理解。
代码
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 10000000
#define X 1000000007
using namespace std;
int n,m,k;
I int Qpow(RI x,RI y) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
template<int SZ> class LinearSieve//线性筛
{
private:
int Pt,P[SZ+5],f[SZ+5],p[SZ+5];
public:
I int operator [] (CI x) Con {return f[x];}
I void Init()
{
RI i,j;for(f[1]=1,i=2;i<=SZ;++i)
{
!P[i]&&(P[++Pt]=i,p[i]=Qpow(i,k),f[i]=p[i]-1);
for(j=1;j<=Pt&&i*P[j]<=SZ;++j)
if(P[i*P[j]]=1,i%P[j]) f[i*P[j]]=1LL*f[i]*f[P[j]]%X;
else {f[i*P[j]]=1LL*f[i]*p[P[j]]%X;break;}
}
for(i=2;i<=SZ;++i) f[i]=(f[i]+f[i-1])%X;
}
};LinearSieve<N> F;
int main()
{
RI Tt,i,j,l,r,t;scanf("%d%d",&Tt,&k),F.Init();W(Tt--)
{
for(scanf("%d%d",&n,&m),t=0,l=1;l<=min(n,m);l=r+1)//除法分块
r=min(n/(n/l),m/(m/l)),t=(1LL*(n/l)*(m/l)%X*(F[r]-F[l-1]+X)+t)%X;
printf("%d
",t);
}return 0;
}