如果这篇博客帮助到你，可以请我喝一杯咖啡~
CC BY-NC-SA 4.0 （除特别声明或转载文章外）
容易发现，这道题在高斯消元完成后
\[m_{i,j}'=\begin{cases}0&i\nmid j\\f(i)&i\mid j\end{cases}\]原因是只有 $i$ 的因数行会影响第 $i$ 行的值，而对于 $i\mid j$ 的 $j$，这些行的对应位置都相等。而对于 $i\nmid j$ 的位置，因为 $\gcd(i,j)\mid i，\gcd(\gcd(i,j),j)=\gcd(i,j)$，所以在消元到第 $\gcd(i,j)$ 行时，$m_{i,j}=m_{\gcd(i,j),j}$，这样一减就变成 $0$ 了。
我们有 $m_{i,i}=\gcd(i,i)^k-\sum\limits_{d\mid i,d\ne i}m_{d,i}$ 即 $f(i)=i^k-\sum\limits_{d\mid i}f(d)$，也就是 $\sum\limits_{d\mid i}f(i)=id_k(i)$，到这一步写个狄利克雷除法就可以 $O(n\log n)$ 的做了，但根据数据范围看，我们应该需要 $O(n)$ 的解法，所以继续推，设 $f$ 的 DGF 为 $F(x)$，我们已知 $1$ 的 DGF 为 $\zeta(x)$，$id_k$ 的 DGF 为 $\zeta(x-k)$，所以 $F(x)\zeta(x)=\zeta(x-k)$，$F(x)=\frac{\zeta(x-k)}{\zeta(x)}$，事实上，$k=1$ 时 $f$ 就是 $\varphi$。
\[\begin{aligned} F(x)=&\frac{\zeta(x-k)}{\zeta(x)}\\ =&\prod\limits_{p\in \mathrm{Prime}} \frac{1-p^{-x}}{1-p^{-x+k}}\\ =&\prod\limits_{p\in \mathrm{Prime}}\left(1+\sum\limits_{i=1}^\infin\frac{p^{ik}-p^{(i-1)k}}{p^{ix}}\right) \end{aligned}\]所以 $f$ 在素数幂 $p^i$ 处的取值为 $p^{ik}-p^{(i-1)k}$，线性筛即可。
答案即为 $\prod\limits_{i=1}^nf(i)$。
#include<cstdio>
int pw[1000010],f[1000010],t,n,k,np[1000010],p[1000010];
int const mod=1e6+3;
int fpow(int x,int y){
int res=1;
while(y){
if(y&1)res=1ll*res*x%mod;
x=1ll*x*x%mod;
y>>=1;
}
return res;
}
int main(){
scanf("%d",&t);
while(t--){
scanf("%d%d",&n,&k);
f[1]=1;
int ans=1,cnt=0;
for(int i=2;i<=n;i++){
if(!np[i])p[++cnt]=i,f[i]=(mod+(pw[i]=fpow(i,k))-1)%mod;
for(int j=1;j<=cnt&&p[j]*i<=n;j++){
np[i*p[j]]=1;
if(i%p[j]) f[i*p[j]]=1ll*f[i]*f[p[j]]%mod;
else{f[i*p[j]]=1ll*f[i]*pw[p[j]]%mod;break;}
}
ans=1ll*ans*f[i]%mod;
}
printf("%d\n",ans);
}
return 0;
}