题面
题解
首先我们将原矩阵写成 \(A+B\)，其中 \(B\) 全是 \(C\)，那么 \(A\) 的第 \(i\) 行就只有其倍数处有值，且 \(A_{i,i}=1-C,A_{i,j(i|j\land i\neq j)}=-C\)。
那么原来的行列式就变成了：
\[\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} \]考虑在 \(A\) 中选的点数（\(|S|\)），发现它不可能小于 \(n-1\)，否则对于 \(B\) 要求的就是一个大小大于 \(1\)、全是 \(C\) 的矩阵的行列式，它恒为 \(0\)。
那么假设 \(B\) 中选的那个点（\(B\) 中没选点的情况很好算，略）为 \((i,p_i)\)。考虑把行列式问题放在图论上（注意到 \((-1)^{\operatorname{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\land d\neq n}\frac{C}{1-C}f_d \]其中前者是链长为 \(1\) 的方案，后者是链长大于 \(1\) 的方案，枚举链尾的前一个数 \(d\)，注意环个数的变化所带来的正负号影响。
我们只需要求出 \(\sum_{i=1}^n f_n\) 即可，发现这也是可以用类似杜教筛的方法来做的：
\[\begin{aligned} \sum_{i=1}^nf_i&=\sum_{i=1}^n\left(C(1-C)^{n-1}\sum_{d|i,d\neq i}\frac{C}{1-C}f_d\right)\\ &=nC(1-C)^{n-1}+\sum_{d=1}^nf_d\cdot\frac{C}{1-C}\left(\left\lfloor\frac{n}{d}\right\rfloor-1\right) \end{aligned} \]有一个问题是如何线性预处理小范围的 \(\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\}\) 选一个符合它的最小的数作为代表元（具体来说，将 \(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
*/