摆
题目描述
题解
首先，我们观察一下这个矩阵，看它有什么性质，你会发现它是长这个样子的：
一个上三角的部分全部都是
C
C
C，中线是
1
1
1，下面有的是
C
C
C有的是
0
0
0。
由于矩阵上把某一行加或减在另一行上，并不会改变该矩阵行列式的值，我们不妨将每一行都减去它下一行的值，于是你会发现它变成这个样子了：
画的好丑呀
也就是一个上海森堡矩阵，右上是一个全
0
0
0的三角，左下角
i
i
i的倍数处为
C
C
C，
i
i
i的倍数减一处为
?
C
-C
?C。
首先对于这种海森堡矩阵，我们的行列式可以考虑从行列式的定义入手求解。
定义的方法相当于是我们枚举一个排列，计算排列的乘积。
由于排列中肯定是存在置换环的，不妨考虑单个置换环会怎么产生贡献。
显然，在这个环中肯定存在一个
i
i
i使得
p
i
>
i
p_i>i
pi?>i，而走到
p
i
p_i
pi?后，我们的
p
p
i
p_{p_i}
ppi??肯定是
p
i
?
1
p_i-1
pi??1了。
因为我们的上边只存在
(
i
?
1
,
i
)
(i-1,i)
(i?1,i)的点，我们往回走就只能从这个点走，再往后走就回不来了，而且往回走还每次只能走一步。
我们考虑这个环会产生怎样的贡献，走回来时是每步贡献
C
?
1
C-1
C?1，走过去时每步贡献
C
C
C，并且每个环还会贡献一个逆序对。
看起来太麻烦了，我们干脆给每个点都除去一个
C
?
1
C-1
C?1，这样就相当于我们一个环贡献就直接乘上一个
C
1
?
C
\frac{C}{1-C}
1?CC?。
我们定义
f
i
f_i
fi?表示大小为
i
×
i
i\times i
i×i的上海森堡矩阵的行列式，容易得到转移方程:
f
i
=
f
i
?
1
+
C
1
?
C
∑
d
∣
i
∧
d
≠
i
f
d
?
f
d
?
1
f_i=f_{i-1}+\frac{C}{1-C}\sum_{d|i\wedge d\neq i}f_d-f_{d-1}
fi?=fi?1?+1?CC?d∣i∧d?=i∑?fd??fd?1?把
f
i
?
1
f_{i-1}
fi?1?减到左边去，就变成了一个差分的形式，记差分的为
g
i
g_i
gi?，有转移式：
g
i
=
∑
d
∣
i
∧
d
≠
i
C
1
?
C
g
d
g_i=\sum_{d|i\wedge d\neq i}\frac{C}{1-C}g_{d}
gi?=d∣i∧d?=i∑?1?CC?gd?显然，
g
g
g的前缀和可以通过杜教筛求解，我们考虑将
C
1
?
C
I
\frac{C}{1-C}I
1?CC?I跟
g
g
g卷在一起，就成了经典的杜教筛形式。
我们可以先预处理出来前
n
2
3
n^{\frac{2}{3}}
n32?处的
g
i
g_{i}
gi?值，后面的部分就可以
O
(
n
2
3
)
O\left(n^{\frac{2}{3}}\right)
O(n32?)的数论分块快速计算。
问题就是前面的
g
g
g值也是不能暴力计算的，否则会
T
T
T飞。
可以考虑对于
m
=
∏
p
i
a
i
m=\prod p_i^{a_i}
m=∏piai??，
g
m
g_m
gm?的值显然只与集合
A
A
A有关。
所以我们不妨考虑对于每种集合，选一个数暴力
O
(
m
)
O\left(\sqrt{m}\right)
O(m?)计算。
显然，这样的集合个数是
O
(
P
(
log
?
m
)
)
O(P(\log m))
O(P(logm))级别的，这部分的计算时间复杂度比较小。
当然，对于集合的存储我们可以考虑
h
a
s
h
hash
hash，这样就能相当快地找到它属于哪个集合了。
加上
H
a
s
h
M
a
p
HashMap
HashMap这部分就能线性了。
不过需要注意的是我们算出来
f
n
f_n
fn?并不一定是答案，还要乘上我们之前除去的
(
1
?
C
)
n
?
1
(1-C)^{n-1}
(1?C)n?1，其中第一行没除，所以是
n
?
1
n-1
n?1次方，我们的答案实际上是
(
1
?
C
)
n
?
1
f
n
(1-C)^{n-1}f_n
(1?C)n?1fn?。
总时间复杂度 O ( n 2 3 ) O\left(n^{\frac{2}{3}}\right) O(n32?)。
源码
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef pair<int,int> pii;
typedef unsigned int uint;
#define MAXN 20000005
#define MAXM 1000010
#define pb push_back
#define mkpr make_pair
#define fir first
#define sec second
const LL INF=0x3f3f3f3f3f3f3f3f;
const int mo=998244353;
const int mod=1e6+7;
const int lim=20000000;
template<typename _T>
void read(_T &x){
_T f=1;x=0;char s=getchar();
while(s<'0'||s>'9'){if(s=='-')f=-1;s=getchar();}
while('0'<=s&&s<='9'){x=(x<<3)+(x<<1)+(s^48);s=getchar();}
x*=f;
}
template<typename _T>
_T Fabs(_T x){return x<0?-x:x;}
LL gcd(LL a,LL b){return !b?a:gcd(b,a%b);}
int add(int x,int y,int p){return x+y<p?x+y:x+y-p;}
void Add(int &x,int y,int p){x=add(x,y,p);}
int qkpow(int a,int s,int p){int t=1;while(s){if(s&1)t=1ll*a*t%p;a=1ll*a*a%p;s>>=1;}return t;}
int C,w,g[MAXN],sum[MAXN],prime[MAXM*2],cntp,F[MAXM],idx;LL n,up;
int minp[MAXN],cnt[MAXN],h[MAXN],hs[30],id1[MAXM],id2[MAXM];
bool oula[MAXN];
struct HashMap{
int val[MAXM],id[MAXM],head[MAXM],nxt[MAXM],tot;
void insert(int ai,int aw){
int pos=ai%mod,now=++tot;id[now]=ai;
nxt[now]=head[pos];head[pos]=now;val[now]=aw;
}
int query(int ai){
int pos=ai%mod,now=head[pos];
while(now&&id[now]!=ai)now=nxt[now];
return val[now];
}
}Mp;
void init(){
for(int i=2;i<=lim;i++){
if(!oula[i])prime[++cntp]=i;
for(int j=1;j<=cntp&&1ll*i*prime[j]<=lim;j++){
oula[i*prime[j]]=1;
minp[i*prime[j]]=j;
if(i%prime[j]==0)break;
}
}
}
int getId(LL x){return x<=n/x?id1[x]:id2[n/x];}
int getF(LL x){
int id=getId(x);if(F[id])return F[id];
for(LL l=2,r;l<=x;l=r+1)
r=x/(x/l),Add(F[id],1ll*(r-l+1)%mo*getF(x/l)%mo,mo);
return F[id]=(1ll*w*F[id]+1)%mo;
}
int main(){
//freopen("bigben.in","r",stdin);
//freopen("bigben.out","w",stdout);
read(n);read(C);if(!C){puts("1");return 0;}
if(C==1){puts(n>2?"0":"1");return 0;}
init();sum[1]=g[1]=h[1]=1;mt19937 e(time(0));up=min(n,(LL)lim);
for(int i=1;i<=25;i++)hs[i]=prime[i];shuffle(hs+1,hs+26,e);
w=1ll*C*qkpow(add(1,mo-C,mo),mo-2,mo)%mo;
for(int i=25;i>1;i--)hs[i]=1ll*hs[i]*qkpow(hs[i-1],mo-2,mo)%mo;
for(int i=2;i<=up;i++){
if(!oula[i])h[i]=hs[cnt[i]=1];
for(int j=1;j<=cntp;j++){
int t=i*prime[j];if(t>up)break;minp[t]=j;
if(i%prime[j]==0){h[t]=1ll*h[i]*hs[cnt[t]=cnt[i]+1]%mo;break;}
else h[t]=1ll*h[i]*hs[cnt[t]=1]%mo;
}
g[i]=Mp.query(h[i]);
if(!g[i]){
g[i]=1;int ni=sqrt(i);
for(int j=2;j<=ni;j++)if(i%j==0){
Add(g[i],g[j],mo);
if(j!=i/j)Add(g[i],g[i/j],mo);
}
g[i]=1ll*w*g[i]%mo;
Mp.insert(h[i],g[i]);
}
sum[i]=add(sum[i-1],g[i],mo);
}
for(LL l=1,r;l<=n;l=r+1){
r=n/(n/l);(n/l<=l?id1[n/l]:id2[l])=++idx;
if(n/l<=lim)F[idx]=sum[n/l];
}
int ans=1ll*qkpow(add(1,mo-C,mo),(n-1)%(mo-1),mo)*getF(n)%mo;
printf("%d\n",ans);
return 0;
}