【NOI模拟赛】摆(线性代数,杜教筛)

题面

6s , 1024mb

我是XYX,我擅长摆。

我在摆大烂的时候看到一个 n × n n\times n n×n 的矩阵 A A A
A i , j = { 1 i = j 0 i ≠ j ∧ i ∣ j C o t h e r w i s e A_{i,j}=\begin{cases} 1 & i=j\\ 0 & i\not=j\land i|j\\ C & {\rm otherwise} \end{cases} Ai,j=10Ci=ji=jijotherwise

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

1 ≤ n ≤ 1 0 11 , 0 ≤ C < 998244353 1\leq n\leq10^{11},0\leq C<998244353 1n1011,0C<998244353

样例

6 2 → \rightarrow 998244350

99 4 → \rightarrow 714389845

10000000000 10 → \rightarrow 82177551

题解

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

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

这是矩阵 A A A
[ 1 0 0 0 0 … C 1 C 0 C … C C 1 C C … C C C 1 C … C C C C 1 … … … … … … ⋱ ] \left[\begin{matrix} 1 & 0 & 0 & 0 &0&\ldots\\ C & 1 & C & 0 & C&\ldots\\ C & C & 1 & C & C&\ldots\\ C & C & C & 1 & C&\ldots\\ C & C & C & C & 1 &\ldots\\ \ldots&\ldots&\ldots&\ldots&\ldots&\ddots \end{matrix}\right] 1CCCC01CCC0C1CC00C1C0CCC1

我们把矩阵 A A A 每一行减去上面一行(从下面开始,每一行减去原来的上面一行,第一行不动),可以得到一个上海森堡矩阵:
[ 1 0 0 0 0 … C − 1 1 C 0 C … 0 C − 1 1 − C C 0 … 0 0 C − 1 1 − C 0 … 0 0 0 C − 1 1 − C … … … … … … ⋱ ] \left[\begin{matrix} 1 & 0 & 0 & 0 &0&\ldots\\ C-1 & 1 & C & 0 & C&\ldots\\ 0 & C-1 & 1-C & C & 0&\ldots\\ 0 & 0 & C-1 & 1-C & 0&\ldots\\ 0 & 0 & 0 & C-1 & 1-C &\ldots\\ \ldots&\ldots&\ldots&\ldots&\ldots&\ddots \end{matrix}\right] 1C100001C1000C1CC1000C1CC10C001C

该矩阵每个置换环都是一段标号连续的逆时针环,类似于 i → j → j − 1 → . . . → i + 1 i\rightarrow j\rightarrow j-1 \rightarrow...\rightarrow i+1 ijj1...i+1 的。差分可得,右上方的每个倍数关系 i ∣ j i|j ij ,会给 A i , j A_{i,j} Ai,j 贡献 − C -C C ,给 A i + 1 , j A_{i+1,j} Ai+1,j 贡献 C C C。我们将矩阵每个位置乘 1 1 − C \frac{1}{1-C} 1C1 ,然后令 f i f_i fi 为保留前 i i i i i i 列的行列式,那么
f i = f i − 1 + C 1 − C ∑ j ∣ i , j < i ( f j − f j − 1 ) f_i=f_{i-1}+\frac{C}{1-C}\sum_{j|i,j<i} (f_j-f_{j-1}) fi=fi1+1CCji,j<i(fjfj1)

g i = f i − f i − 1 g_i=f_i-f_{i-1} gi=fifi1 ,那么 g i = C 1 − C ∑ j ∣ i , j < i g j g_i=\frac{C}{1-C}\sum_{j|i,j<i}g_j gi=1CCji,j<igj

我们要求的是 g i g_i gi 的前缀和 f n f_n fn ,看数据范围应该是个亚线性筛。

h 1 = − 1 , h i = C 1 − C h_1=-1,h_{i}=\frac{C}{1-C} h1=1,hi=1CC ,那么(狄利克雷卷积) h ∗ g = − e h*g=-e hg=e ,所以可以杜教筛(杜教筛并没限制只能是积性函数)。
∑ i = 1 n ∑ j ∣ i h j g i / j = ∑ j n h j ∑ i n / j g i = − 1 ⇒ f n = 1 + ∑ j = 2 n h j f n / j \sum_{i=1}^{n}\sum_{j|i}h_jg_{i/j}=\sum_{j}^nh_j\sum_{i}^{n/j}g_i=-1\\ \Rightarrow f_n=1+\sum_{j=2}^nh_jf_{n/j} i=1njihjgi/j=jnhjin/jgi=1fn=1+j=2nhjfn/j

我们预处理 n 2 / 3 n^{2/3} n2/3 以内的 f i f_i fi ,那么杜教筛的复杂度就是 O ( n 2 / 3 ) O(n^{2/3}) O(n2/3)

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

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

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

总复杂度 O ( n 2 / 3 ) O(n^{2/3}) O(n2/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;
}
  • 2
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值