【SP1772 DETER2】Find The Determinant II （线性代数，狄利克雷逆）

题面

对于一个 $N\times N$ 的矩阵，$m_{i,j}=\gcd (i,j{)}^k$

你的任务是计算这个矩阵的行列式

共 $t\le 20$ 组数据，$1\le N\le 10^6,1\le k\le 10^9$ 。

题解

看这个数据范围，应该要我们达到预处理不超过 $O(N\log N)$ ，循环内不超过 $O(N)$ 的复杂度。

如果做过 Find The Determinant（$m_{i,j}=gcd(i,j)$），应该知道，我们其实是利用了 $n={\sum }_{i|n}\varphi (i)$ ，然后将 $M$ 拆成两个矩阵相乘：$m_{i,j}=gcd(i,j)={\sum }_{k|i,k|j}\varphi (k)={\sum }_k([k|i]\varphi (k))\cdot [k|j]$ ，利用这两个矩阵都是三角矩阵从而直接累乘，得出最终行列式答案：$\prod \varphi (i)$。

这道题也是一样的思路，不过 $n^k$ 又等于什么函数的狄利克雷卷积呢？

我们假设已经知道了这个函数，令 $n^k={\sum }_{i|n}f(i)$ ，那么
$$m_{i,j}=gcd(i,j{)}^k=\sum _{k|i,k|j}f(k)=\sum _k([k|i]f(k))\cdot [k|j]$$

我们就可以得到，最终的答案是 ${\prod }_{i=1}^{N}f(i)$ 。

问题是现在不知道 $f(i)$ 每个位置的函数值，我们需要求。

$n^k={\sum }_{i|n}f(i)$ 其实可以写成狄利克雷卷积的形式：$i{d}_k=f\ast I$ 。

也就是说，$f=i{d}_k\ast I^{-1}$ ，我们可以求狄利克雷逆！

别急，我们只需要知道有狄利克雷逆这个东西就行了，然后可以现推：
$$\begin{gathered} I^{-1}\ast I=e \\ {I}^{-1}(1)I(1)=1\Rightarrow I^{-1}(1)=1 \\ (n>1)\sum _{i|n}{I}^{-1}(i)I(n/i)=0 \\ \Rightarrow I^{-1}(n)I(1)+\sum _{i|n,i<n}{I}^{-1}(i)I(n/i)=0 \\ \Rightarrow I^{-1}(n)=-\sum _{i|n,i<n}{I}^{-1}(i)I(n/i) \end{gathered}$$

可以用埃氏筛法模拟这个递推过程，时间复杂度 $O(N\log N)$ ，由于 $I^{-1}$ 并不随着 $N$ 和 $k$ 变化，因此我们可以在所有询问之前预处理出来 $I^{-1}(1\sim 10^6)$。

---

接下来我们需要线性时间复杂度求出每个 $f(i)={\sum }_{j|i}i{d}_k(j)I^{-1}(i/j)$ 。

这里有个非常好理解但是证明稍麻烦的结论：积性函数 $g$ 的狄利克雷逆 $g^{-1}$ 也是积性函数。

所以，$I^{-1}$ 是积性函数，进而 $f$ 也是积性函数。$f$ 是积性函数，我们便可以思考欧拉筛。

由于 $i{d}_k$ 完全积性，我们对于每个质数 $p$ 用快速幂求出 $i{d}_k(p)$，再欧拉筛可以求出所有 $i{d}_k$ 。再对于每个 $p^w$（质数的幂）枚举因数（$1,p,p^2,...,p^w$，无需根号）求出 $f(p^w)$ ，进而欧拉筛求出所有 $f$ 。时间复杂度 $O(\frac{N\log k}{\ln N})\approx O(N)$ 。

具体实现上，我们可以欧筛的时候求每个数最小质因子的乘积，用于判断是否是质数幂，以及用于求积性函数。

CODE

注意，这个屑OJ限制了代码长度

#include<map>
#include<set>
#include<cmath>
#include<queue>
#include<stack>
#include<random>
#include<vector>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define MAXN 1000005
#define LL long long
#define ULL unsigned long long
#define ENDL putchar('\n')
#define DB long double
#define lowbit(x) (-(x) & (x))
#define FI first
#define SE second
LL read() {
	LL f = 1,x = 0;int s = getchar();
	while(s < '0' || s > '9') {if(s<0)return -1;if(s=='-')f=-f;s = getchar();}
	while(s >= '0' && s <= '9') {x = (x<<1) + (x<<3) + (s^48);s = getchar();}
	return f*x;
}
const int MOD = 1000003;
int n,m,s,o,k;
int iv[MAXN];
int p[MAXN],f[MAXN],cnt,f1[MAXN],ls[MAXN];
int fn[MAXN];
int qkpow(int a,int b) {
	int res = 1; while(b > 0) {
		if(b & 1) res = res *1ll* a % MOD;
		a = a *1ll* a % MOD; b >>= 1;
	}return res;
}
void sieve(int n) {
	f1[1] = 1; fn[1] = 1; cnt = 0;
	for(int i = 1;i <= n;i ++) f[i] = 0;
	for(int i = 2;i <= n;i ++) {
		if(!f[i]) p[++ cnt] = i,f1[i] = qkpow(i,k),fn[i] = (f1[i] + iv[i]) % MOD,ls[i] = i;
		for(int j = 1;j <= cnt && i*p[j] <= n;j ++) {
			int y = i*p[j];
			f[y] = 1; 
			f1[y] = f1[i]*1ll*f1[p[j]] % MOD;
			if(i % p[j] == 0) {
				ls[y] = ls[i]*p[j];
				fn[y] = 0;
				if(y == ls[y]) {
					for(int pw = 1;pw <= ls[y];pw *= p[j]) (fn[y] += f1[pw]*1ll*iv[y/pw]%MOD) %= MOD;
				}
				else fn[y] = fn[y/ls[y]]*1ll*fn[ls[y]] % MOD;
				break;
			}
			ls[y] = p[j];
			fn[y] = fn[i]*1ll*fn[p[j]] % MOD;
		}
	}return ;
}
int main() {
	iv[1] = 1;
	for(int i = 1;i <= MAXN-5;i ++) {
		for(int j = 2;i*j <= MAXN-5;j ++) {
			(iv[i*j] += MOD - iv[i]) %= MOD;
		}
	}
	int T = read();
	while(T --) {
		n = read(); k = read();
		sieve(n);
		int ans = 1;
		for(int i = 1;i <= n;i ++) {
			ans = ans *1ll* fn[i] % MOD;
		}
		printf("%lld\n",ans);
	}
	return 0;
}