题意
(T) 组数据，每组给定两个整数 (n,k)，求 (det A)，其中 (A) 为一个 (n imes n) 的矩阵且 (A_{i,j}=gcd(i,j)^k)，对 (10^6+3) 取模。
( exttt{Data Range:}1leq Tleq 20,1leq nleq 10^6,1leq kleq 10^9)
题解
很好的一道数论题。（然而可能只是因为我出过一道相似的）
类似于这个题的思路，设 (f_{i}) 为消元后 (A_{i,i}) 的值，我们有
[f_n=n^k-sumlimits_{dmid n,d<n}f_d
]
移项得到
[n^k=sumlimits_{dmid n}f_d
]
这个式子跟 (varphi(n)) 的狄利克雷卷积式有些类似，我们考虑狄利克雷生成函数求通项。为了不失普遍性，我们记消元之后的 (f(n)=varphi_k(n))。（这个东西其实叫 Jordan totient function，符号为 (J_k(n))，但是为了体现出类比所以我选用这个符号）
套个 (sum) 我们有
[sumlimits_{igeq 1}frac{i^k}{i^x}=zeta(x)sumlimits_{igeq 1}frac{varphi_k(i)}{i^x}
]
移项有
[sumlimits_{igeq 1}frac{varphi_k(i)}{i^x}=frac{zeta(x-k)}{zeta(x)}
]
使用欧拉乘积公式化开右边（其中 (P) 是素数集）
[sumlimits_{igeq 1}frac{varphi_k(i)}{i^x}=prod_{pin P}frac{1-p^{-x}}{1-p^{k-x}}
]
这里有一个定理：如果 (g_n) 存在积性，那么可以将 (g_n) 对应的狄利克雷生成函数写成如下形式
[G(x)=prod_{pin P}left(sum_{i}frac{g_{p^i}}{p^{ix}}
ight)
]
发现之前式子的右边很像这个东西，所以右边写成这个形式有
[sumlimits_{igeq 1}frac{varphi_k(i)}{i^x}=prod_{pin P}left(sum_{i}frac{p^{ik}-p^{(i-1)k}}{p^{ix}}
ight)
]
所以我们知道 (varphi_k(i)) 是积性函数，也知道在 (p^k) 处的取值，所以有：
[varphi_k(n)=n^kprod_{i=1}^{m}frac{p_i^k-1}{p_i^k}
]
线性筛出这东西即可，答案为 (prodlimits_{i=1}^{n}varphi_k(i))，时间复杂度 (O(Tn))。
交 SPOJ 的题是可以正常开 pragma 的
代码
#include<bits/stdc++.h>
#pragma GCC optimize("Ofast,unroll-loops")
using namespace std;
typedef int ll;
typedef long long int li;
const ll MAXN=1e6+51,MOD=1e6+3;
ll test,n,kk,ptot,res;
ll np[MAXN],prime[MAXN],phi[MAXN],pwPr[MAXN];
inline ll read()
{
register ll num=0,neg=1;
register char ch=getchar();
while(!isdigit(ch)&&ch!='-')
{
ch=getchar();
}
if(ch=='-')
{
neg=-1;
ch=getchar();
}
while(isdigit(ch))
{
num=(num<<3)+(num<<1)+(ch-'0');
ch=getchar();
}
return num*neg;
}
inline ll qpow(ll base,ll exponent)
{
ll res=1;
while(exponent)
{
if(exponent&1)
{
res=(li)res*base%MOD;
}
base=(li)base*base%MOD,exponent>>=1;
}
return res;
}
inline void sieve(ll limit)
{
np[1]=phi[1]=1,ptot=0;
for(register int i=2;i<=limit;i++)
{
if(!np[i])
{
prime[++ptot]=i,phi[i]=(pwPr[i]=qpow(i,kk))-1;
}
for(register int j=1;i*prime[j]<=limit;j++)
{
np[i*prime[j]]=1;
if(i%prime[j]==0)
{
phi[i*prime[j]]=(li)phi[i]*pwPr[prime[j]]%MOD;
break;
}
phi[i*prime[j]]=(li)phi[i]*phi[prime[j]]%MOD;
}
}
}
inline void solve()
{
n=read(),kk=read()%(MOD-1),sieve(n),res=1;
for(register int i=1;i<=n;i++)
{
res=(li)res*phi[i]%MOD;
}
printf("%d
",res);
}
int main()
{
test=read();
for(register int i=0;i<test;i++)
{
solve();
}
}