【NOI模拟赛】摆（线性代数，杜教筛）

题面

6s , 1024mb

我是XYX，我擅长摆。

我在摆大烂的时候看到一个 $n\times n$ 的矩阵 $A$ ：
$$A_{i,j}=\{\begin{array}{ll}1& i=j\\ 0& i\ne j\wedge i|j\\ C& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}$$

现在我想知道 $A$ 的行列式，由于我摆了，所以答案只需要对 $998244353$ 取模。

$1\le n\le 10^{11},0\le C<998244353$ 。

样例：

6 2$\to$998244350

99 4$\to$714389845

10000000000 10$\to$82177551

题解

这道题告诉我们，不一定非得构造上三角或下三角矩阵再算行列式。

我们还可以相对容易地构造海森堡矩阵。

这是矩阵 $A$ ：
$$\left[\begin{array}{cccccc}1& 0& 0& 0& 0& \dots \\ C& 1& C& 0& C& \dots \\ C& C& 1& C& C& \dots \\ C& C& C& 1& C& \dots \\ C& C& C& C& 1& \dots \\ \dots & \dots & \dots & \dots & \dots & \ddots \end{array}\right]$$

我们把矩阵 $A$ 每一行减去上面一行（从下面开始，每一行减去原来的上面一行，第一行不动），可以得到一个上海森堡矩阵：
$$\left[\begin{array}{cccccc}1& 0& 0& 0& 0& \dots \\ C-1& 1& C& 0& C& \dots \\ 0& C-1& 1-C& C& 0& \dots \\ 0& 0& C-1& 1-C& 0& \dots \\ 0& 0& 0& C-1& 1-C& \dots \\ \dots & \dots & \dots & \dots & \dots & \ddots \end{array}\right]$$

该矩阵每个置换环都是一段标号连续的逆时针环，类似于 $i\to j\to j-1\to ...\to i+1$ 的。差分可得，右上方的每个倍数关系 $i|j$ ，会给 $A_{i,j}$ 贡献 $-C$ ，给 $A_{i+1,j}$ 贡献 $C$。我们将矩阵每个位置乘 $\frac{1}{1-C}$ ，然后令 $f_i$ 为保留前 $i$ 行 $i$ 列的行列式，那么
$$f_i=f_{i-1}+\frac{C}{1-C}\sum _{j|i,j<i}(f_j-f_{j-1})$$

设 $g_i=f_i-f_{i-1}$ ，那么 $g_i=\frac{C}{1-C}{\sum }_{j|i,j<i}{g}_j$ 。

我们要求的是 $g_i$ 的前缀和 $f_n$ ，看数据范围应该是个亚线性筛。

令 $h_1=-1,h_i=\frac{C}{1-C}$ ，那么（狄利克雷卷积） $h\ast g=-e$ ，所以可以杜教筛（杜教筛并没限制只能是积性函数）。
$$\begin{gathered} \sum _{i=1}^n\sum _{j|i}{h}_j{g}_{i/j}=\sum _j^n{h}_j\sum _i^{n/j}{g}_i=-1 \\ \Rightarrow f_n=1+\sum _{j=2}^n{h}_j{f}_{n/j} \end{gathered}$$

我们预处理 $n^{2/3}$ 以内的 $f_i$ ，那么杜教筛的复杂度就是 $O(n^{2/3})$ 。

但是预处理很容易带个 $\log$ ，所以需要优化。

令 $x=\prod p_i^{w_i}$ ，我们把 $p_i$ 从小到大依次换成 $2,3,5,7,...$ 得到 $y=\prod p{{}^{\prime }}_i^{w_i}$ ，容易发现 $g_x=g_y$ ，因为 $g$ 只跟 $w_i$ 的可重集有关。我们用欧拉筛求出每个数的最大质因子和质因子种数，进而递推求出每个 $x$ 对应的 $y$ 。暴力发现本质不同的 $y$ 只有最多一千多个，我们可以每个 $O(n^{1/3})$ 求出 $g$ 。

注意，由于矩阵每个位置乘了 $\frac{1}{1-C}$ ，且由于矩阵左上角并不标准，所以最终答案是
$$f_n\cdot (1-C{)}^{n-1}$$

总复杂度 $O(n^{2/3})$ 。

CODE

#include<map>
#include<set>
#include<cmath>
#include<ctime>
#include<queue>
#include<stack>
#include<random>
#include<bitset>
#include<vector>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<unordered_map>
#pragma GCC optimize(2)
using namespace std;
#define MAXN 10000005
#define LL long long
#define ULL unsigned long long
#define ENDL putchar('\n')
#define DB double
#define lowbit(x) (-(x) & (x))
#define FI first
#define SE second
#define PR pair<int,int>
int xchar() {
    
	static const int maxn = 1000000;
	static char b[maxn];
	static int pos = 0,len = 0;
	if(pos == len) pos = 0,len = fread(b,1,maxn,stdin);
	if(pos == len) return -1;
	return b[pos ++];
}
// #define getchar() xchar()
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;
}
void putpos(LL x) {
    if(!x)return ;putpos(x/10);putchar((x%10)^48);}
void putnum(LL x) {
    
	if(!x) {
    putchar('0');return ;}
	if(x<0) putchar('-'),x = -x;
	return putpos(x);
}
void AIput(LL x,int c) {
    putnum(x);putchar(c);}

const int MOD = 998244353;
int n,m,s,o,k;
LL N;
inline 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;
}
int V;
int sg[MAXN],lw[MAXN],ct[MAXN],hb[MAXN];
int p[MAXN],cnt; bool f[MAXN];
void sieve(int n) {
    
	sg[1] = 1; lw[1] = 1; int nm = 0;
	for(int i = 2;i <= n;i ++) {
    
		if(!f[i]) p[++ cnt] = i,lw[i] = 2,hb[i] = i,ct[i] = 1;
		if(lw[i] == i) {
    
			for(int j = 1;j * j <= i;j ++) {
    
				if(i % j == 0) {
    
					(sg[i] += sg[j]) %= MOD;
					if(i/j > j && j > 1) {
    
						(sg[i] += sg[i/j]) %= MOD;
					}
				}
			}
			sg[i] = sg[i] *1ll* V % MOD;
		}
		else sg[i] = sg[lw[i]];
		for(int j = 1;j <= cnt && (nm=i*p[j]) <= n;j ++) {
    
			f[nm] = 1; hb[nm] = hb[i];
			if(i % p[j] == 0) {
    
				ct[nm] = ct[i];
				lw[nm] = lw[nm/hb[nm]] * p[ct[nm]];
				break;
			}
			ct[nm] = ct[i] + 1;
			lw[nm] = lw[nm/hb[nm]] * p[ct[nm]];
		}
	}
	for(int i = 2;i <= n;i ++) {
    
		(sg[i] += sg[i-1]) %= MOD;
	}
	return ;
}
int sh(LL x) {
    
	if(x<1) return 0;
	return ((x-1)%MOD*V%MOD +MOD-1) % MOD;
}
int s1[1000005];
int Sg(LL x) {
    
	if(x <= 10000000) return sg[x];
	int &rs = s1[N/x];
	if(rs >= 0) return rs;
	rs = 1;
	for(LL i = 2;i <= x;i ++) {
    
		LL r = x/(x/i);
		rs += (r-i+1)%MOD*V%MOD *1ll* Sg(x/i) % MOD;
		if(rs >= MOD) rs -= MOD;
		i = r;
	} return rs;
}
int main() {
    
	freopen("bigben.in","r",stdin);
	freopen("bigben.out","w",stdout);
	N = read(); int C = read();
	if(C <= 1) return AIput(N<=2,'\n'),0;
	V = C *1ll* qkpow(MOD+1-C,MOD-2) % MOD;
	sieve(10000000);
	memset(s1,-1,sizeof(s1));
	int as = Sg(N) *1ll* qkpow(MOD+1-C,(N-1)%(MOD-1)) % MOD;
	AIput(as,'\n');
	return 0;
}