【XSY4350】摆（行列式，数论，杜教筛）

题面

摆

题解

首先我们将原矩阵写成 $A+B$，其中 $B$ 全是 $C$，那么 $A$ 的第 $i$ 行就只有其倍数处有值，且 $A_{i,i}=1-C,A_{i,j(i|j\wedge i\ne j)}=-C$。

那么原来的行列式就变成了：
∑ p ( − 1 ) sgn ⁡ ( p ) ∑ S ⊆ { 1 , ⋯   , n } ∏ i ∈ S A i , p i ∏ i ∉ S B i , p i \sum_{p}(-1)^{\operatorname{sgn}(p)}\sum_{S\subseteq \{1,\cdots,n\}}\prod_{i\in S}A_{i,p_i}\prod_{i\not\in S}B_{i,p_i} p∑​(−1)sgn(p)S⊆{1,⋯,n}∑​i∈S∏​Ai,pi​​i​∈S∏​Bi,pi​​
考虑在 $A$ 中选的点数（$|S|$），发现它不可能小于 $n-1$，否则对于 $B$ 要求的就是一个大小大于 $1$、全是 $C$ 的矩阵的行列式，它恒为 $0$。

那么假设 $B$ 中选的那个点（$B$ 中没选点的情况很好算，略）为 $(i,p_i)$。考虑把行列式问题放在图论上（注意到 $(-1{)}^{\mathrm{sgn}(p)}$ 为 $(-1{)}^{n-\text{环个数}}$，所以系数的问题也是可以处理的），相当于说 $B$ 在图上钦定了一条边 $i\to p_i$ 必须选，那么在 $A$ 中选的就应该是：

• $p_i\to \cdots \to i$ 的一条链，其中每一条边后者是前者的倍数，且边权为 $-C$。
• 发现剩下的点在 $A$ 中都只能选自环（从大往小推每个点即可），边权为 $1-C$。

假设 $f_n$ 表示链尾是 $n$ 的所有方案的和，可以得到：
$$f_n=C(1-C{)}^{n-1}+\sum _{d|n\wedge d\ne n}\frac{C}{1-C}{f}_d$$
其中前者是链长为 $1$ 的方案，后者是链长大于 $1$ 的方案，枚举链尾的前一个数 $d$，注意环个数的变化所带来的正负号影响。

我们只需要求出 ${\sum }_{i=1}^n{f}_n$ 即可，发现这也是可以用类似杜教筛的方法来做的：
$$\begin{array}{rl}\sum _{i=1}^n{f}_i& =\sum _{i=1}^n\left(C(1-C{)}^{n-1}\sum _{d|i,d\ne i}\frac{C}{1-C}{f}_d\right)\\ & =nC(1-C{)}^{n-1}+\sum _{d=1}^n{f}_d\cdot \frac{C}{1-C}\left(\lfloor \frac{n}{d}\rfloor -1\right)\end{array}$$
有一个问题是如何线性预处理小范围的 ${\sum }_{i=1}^n{f}_i$。发现 $f_i$ 只与 $i$ 的因子的形态有关：具体来说，设 $i=p_1^{a_1}\cdots p_k^{a_k}$，那么 $f_i$ 只与可重集 { a i } \{a_i\} {ai​} 有关。

那么我们为每种 { a i } \{a_i\} {ai​} 选一个符合它的最小的数作为代表元（具体来说，将 $a_i$ 从大到小排序后，取代表元为 $2^{a_1}{3}^{a_2}\cdots$），并只对这些代表元暴力求。

代表元数量很小（实测小于等于 $2\times 10^7$ 的代表元只有不到 $600$ 个，所以对它们每个都根号暴力求都没问题），这样就能做到线性预处理 $O(n^{\frac{2}{3}})$ 内的前缀和了。

时间复杂度 $O(n^{\frac{2}{3}})$。

#include<bits/stdc++.h>

#define ll long long

using namespace std;

namespace modular
{
	const int mod=998244353;
	inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
	inline int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
	inline int mul(int x,int y){return 1ll*x*y%mod;}
	inline void Add(int &x,int y){x=x+y>=mod?x+y-mod:x+y;}
	inline void Dec(int &x,int y){x=x-y<0?x-y+mod:x-y;}
	inline void Mul(int &x,int y){x=1ll*x*y%mod;}
}using namespace modular;

inline int poww(int a,ll b)
{
	int ans=1;
	while(b)
	{
		if(b&1) ans=mul(ans,a);
		a=mul(a,a);
		b>>=1;
	}
	return ans;
}

inline int read()
{
	int x=0,f=1;
	char ch=getchar();
	while(ch<'0'||ch>'9')
	{
		if(ch=='-') f=-1;
		ch=getchar();
	}
	while(ch>='0'&&ch<='9')
	{
		x=(x<<1)+(x<<3)+(ch^'0');
		ch=getchar();
	}
	return x*f;
}

const int n2=300000,n3=20000000;
const int N2=n2+10,N3=n3+10;

ll n;
int C,K,coef;
int idx,id1[N2],id2[N2];
int cnt,prime[N3],mup[N3],omega[N3],sign[N3];
bool notprime[N3];
int sumf[N2],pref[N3],ff[N3];

int &id(ll x){return x<=n2?id1[x]:id2[n/x];}

void init()
{
	mup[1]=1,omega[1]=0;
	for(int i=2;i<=n3;i++)
	{
		if(!notprime[i])
		{
			prime[++cnt]=i;
			mup[i]=i,omega[i]=1;
		}
		for(int j=1,v;j<=cnt&&(v=i*prime[j])<=n3;j++)
		{
			notprime[v]=1;
			if(!(i%prime[j]))
			{
				mup[v]=mup[i],omega[v]=omega[i];
				break;
			}
			mup[v]=mup[i]*prime[j],omega[v]=omega[i]+1;
		}
	}
	vector<int> sp;
	sp.push_back(1);
	for(int i=1;i<=cnt;i++)
	{
		sp.push_back(sp.back()*prime[i]);
		if(sp.back()>n3) break;
	}
	sign[1]=1,ff[1]=K;
	for(int i=2;i<=n3;i++)
	{
		sign[i]=sign[i/mup[i]]*sp[omega[i]];
		if(i!=sign[i])
		{
			ff[i]=ff[sign[i]];
			continue;
		}
		ff[i]=K;
		for(int j=1;j*j<=i;j++)
		{
			if(!(i%j))
			{
				Add(ff[i],mul(coef,ff[j]));
				if(j>1&&j*j!=i) Add(ff[i],mul(coef,ff[i/j]));
			}
		}
	} 
	for(int i=1;i<=n3;i++) pref[i]=add(pref[i-1],ff[i]);
}

int solve(ll n)
{
	if(n<=n3) return pref[n];
	if(sumf[id(n)]!=-1) return sumf[id(n)];
	int ans=mul(n%mod,K);
	for(ll l=1,r;;l=r+1)
	{
		r=n/(n/l);
		if(r==n) break;
		Add(ans,mul(mul(coef,(n/l-1)%mod),dec(solve(r),solve(l-1))));
	}
	return sumf[id(n)]=ans;
}

int main()
{
	cin>>n>>C;
	if(C==1)
	{
		puts(n<=2?"1":"0");
		return 0;
	}
	K=mul(C,poww(dec(1,C),n-1));
	coef=mul(C,poww(dec(1,C),mod-2));
	for(ll l=1,r;l<=n;l=r+1)
		r=n/(n/l),id(n/l)=++idx;
	init();
	memset(sumf,-1,sizeof(sumf));
	int res=solve(n);
	printf("%d\n",add(poww(dec(1,C),n),res));
	return 0;
}
/*
2 2
*/