SP1772 Find The Determinant II
阿新 • • 發佈:2020-10-07
題意
\(T\) 組資料,每組給定兩個整數 \(n,k\),求 \(\det A\),其中 \(A\) 為一個 \(n\times n\) 的矩陣且 \(A_{i,j}=\gcd(i,j)^k\),對 \(10^6+3\) 取模。
\(\texttt{Data Range:}1\leq T\leq 20,1\leq n\leq 10^6,1\leq k\leq 10^9\)
題解
很好的一道數論題。(然而可能只是因為我出過一道相似的)
類似於這個題的思路,設 \(f_{i}\) 為消元后 \(A_{i,i}\) 的值,我們有
\[f_n=n^k-\sum\limits_{d\mid n,d<n}f_d \]移項得到
\[n^k=\sum\limits_{d\mid n}f_d \]這個式子跟 \(\varphi(n)\) 的狄利克雷卷積式有些類似,我們考慮狄利克雷生成函式求通項。為了不失普遍性,我們記消元之後的 \(f(n)=\varphi_k(n)\)。(這個東西其實叫 Jordan totient function,符號為 \(J_k(n)\),但是為了體現出類比所以我選用這個符號)
套個 \(\sum\) 我們有
\[\sum\limits_{i\geq 1}\frac{i^k}{i^x}=\zeta(x)\sum\limits_{i\geq 1}\frac{\varphi_k(i)}{i^x} \]移項有
\[\sum\limits_{i\geq 1}\frac{\varphi_k(i)}{i^x}=\frac{\zeta(x-k)}{\zeta(x)} \]使用尤拉乘積公式化開右邊(其中 \(P\) 是素數集)
\[\sum\limits_{i\geq 1}\frac{\varphi_k(i)}{i^x}=\prod_{p\in P}\frac{1-p^{-x}}{1-p^{k-x}} \]這裡有一個定理:如果 \(g_n\) 存在積性,那麼可以將 \(g_n\) 對應的狄利克雷生成函式寫成如下形式
\[G(x)=\prod_{p\in P}\left(\sum_{i}\frac{g_{p^i}}{p^{ix}}\right) \]發現之前式子的右邊很像這個東西,所以右邊寫成這個形式有
\[\sum\limits_{i\geq 1}\frac{\varphi_k(i)}{i^x}=\prod_{p\in P}\left(\sum_{i}\frac{p^{ik}-p^{(i-1)k}}{p^{ix}}\right) \]所以我們知道 \(\varphi_k(i)\) 是積性函式,也知道在 \(p^k\) 處的取值,所以有:
\[\varphi_k(n)=n^k\prod_{i=1}^{m}\frac{p_i^k-1}{p_i^k} \]線性篩出這東西即可,答案為 \(\prod\limits_{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\n",res);
}
int main()
{
test=read();
for(register int i=0;i<test;i++)
{
solve();
}
}