好题啊!可惜没有链接!
题目描述
题解
计算方法和正解是一样的,但是推导过程用不着海森堡矩阵。
考虑把对角线上的元素全部替换为多项式
C
+
x
C+x
C+x,那么此矩阵的行列式就会变为一个多项式,可以代入
x
=
1
−
C
x=1-C
x=1−C 求出答案。
组合一下选取的
x
x
x 可以得到:
det
(
A
)
=
∑
S
⊂
{
1
,
2
,
.
.
.
,
n
}
det
(
B
S
)
x
n
−
∣
S
∣
\det(A)=\sum_{S\subset\{1,2,...,n\}}\det(B_S)x^{n-|S|}
det(A)=S⊂{1,2,...,n}∑det(BS)xn−∣S∣其中
B
S
B_S
BS 是把
A
A
A 的对角线换成
C
C
C 然后限制在
S
S
S 这个子集上的矩阵。比如:
观察
B
S
B_S
BS 可以发现,如果
S
S
S 中的元素可以两两整除(也就是排序后,前一个整除后一个),那么
B
S
B_S
BS 是一个下三角矩阵,行列式为
C
∣
S
∣
C^{|S|}
C∣S∣;
否则如果不能两两整除,那么一定有
det
(
B
S
)
=
0
\det(B_S)=0
det(BS)=0。这个可以递归证明:
如果
S
S
S 中存在两个数不为其它任何数的因子,那么两个数所在的行都是全
C
C
C,那么矩阵一定不满秩;
反之,
S
S
S 中最大的数一定为所有数的倍数,那么它所在的列只有最后一行是
C
C
C,于是可以往删去最后一行一列的子矩阵递归下去。容易发现最后一定会出现上面不满秩的情况。
我们令
r
=
C
1
−
C
r=\frac{C}{1-C}
r=1−CC,那么问题就转化为求满足元素两两整除的集合
S
S
S 的贡献和:
A
n
s
=
(
1
−
C
)
n
∑
S
元
素
两
两
整
除
r
∣
S
∣
Ans=(1-C)^n\sum_{S元素两两整除}r^{|S|}
Ans=(1−C)nS元素两两整除∑r∣S∣
令
g
x
=
∑
max
{
i
∣
i
∈
S
}
=
x
r
∣
S
∣
g_x=\sum_{\max\{i|i\in S\}=x}r^{|S|}
gx=∑max{i∣i∈S}=xr∣S∣(特别地,
g
1
g_1
g1 表示
S
=
{
1
}
S=\{1\}
S={1} 和空集两种情况),那么可以枚举去掉最大值的集合,得到转移式子
g
i
=
r
∑
j
∣
i
g
j
g_i=r\sum_{j|i}g_j
gi=rj∣i∑gj以及初值
g
1
=
r
+
1
g_1=r+1
g1=r+1。
我们的目的是求出
g
x
g_x
gx 的前缀和,而这个只是看起来像狄利克雷卷积的式子不能用来 Min_25 筛,需要另找他法。
令
s
(
n
)
=
∑
i
=
1
n
g
i
s(n)=\sum_{i=1}^ng_i
s(n)=∑i=1ngi,我们其实可以直接枚举
S
S
S 中的最大值和次大值的倍数关系得到状态转移:
s
(
n
)
=
1
+
r
+
∑
i
=
2
n
s
(
⌊
n
i
⌋
)
s(n)=1+r+\sum_{i=2}^ns(\lfloor\frac{n}{i}\rfloor)
s(n)=1+r+i=2∑ns(⌊in⌋)这样我们就得到了一个直接根号分治枚举的做法。
然而这个做法复杂度为 O ( n 3 4 ) O(n^{\frac{3}{4}}) O(n43),只能过80分,还要优化。这里我们可以沿用杜教筛中的思想,先预处理出 g i g_i gi 的前 n 2 3 n^{\frac{2}{3}} n32 个值求出前缀和,然后剩下的部分根号分治枚举,这样就可以把根号分治的复杂度降至 O ( n 2 3 ) O(n^{\frac{2}{3}}) O(n32)。
现在问题在于,怎么 O ( n ) O(n) O(n) 线性地求出 g i g_i gi。考虑 x x x 的唯一分解,我们发现若所有质因子的次数构成的可重集相同,那么 g x g_x gx 相等(因为这本质上是一个往倍数走的路径计数,所以所有质因子等价)。我们线性筛求出所有元素的可重集的哈希值,那么哈希值相等的就只用求一次。
神奇的是,可重集的种类数非常少,少于 线 性 筛 范 围 \sqrt{线性筛范围} 线性筛范围,所以我们挑同一哈希值中最小的数直接枚举因子暴力求。
最后总复杂度即为 O ( n 2 3 ) O(n^{\frac{2}{3}}) O(n32)。
代码
#include<bits/stdc++.h>//JZM yyds!!
#define ll long long
#define lll __int128
#define uns unsigned
#define fi first
#define se second
#define IF (it->fi)
#define IS (it->se)
#define END putchar('\n')
#define lowbit(x) ((x)&-(x))
#define inline jzmyyds
using namespace std;
const int MAXN=2e6+5;
const ll INF=1e18;
ll read(){
ll x=0;bool f=1;char s=getchar();
while((s<'0'||s>'9')&&s>0){if(s=='-')f^=1;s=getchar();}
while(s>='0'&&s<='9')x=(x<<1)+(x<<3)+(s^48),s=getchar();
return f?x:-x;
}
int ptf[50],lpt;
void print(ll x,char c='\n'){
if(x<0)putchar('-'),x=-x;
ptf[lpt=1]=x%10;
while(x>9)x/=10,ptf[++lpt]=x%10;
while(lpt>0)putchar(ptf[lpt--]^48);
if(c>0)putchar(c);
}
const ll MOD=998244353;
ll ksm(ll a,ll b,ll mo){
ll res=1;
for(;b;b>>=1,a=a*a%mo)if(b&1)res=res*a%mo;
return res;
}
unordered_map<ll,ll>mp;
mt19937 Rand(114514);
ll n,B,C,r,f[MAXN];
ll g[MAXN*10],mi[MAXN*10],hs[MAXN*10],ad[233];
bool nop[MAXN*10];
int pr[MAXN],le,m;
void sieve(int n){
nop[0]=nop[1]=1;
for(int a=2;a<=n;a++){
if(!nop[a])pr[++le]=a,mi[a]=1,hs[a]=ad[1];
for(int i=1,u;i<=le&&(u=pr[i]*a)<=n;i++){
nop[u]=1;
if(a%pr[i]==0){
mi[u]=mi[a]+1,hs[u]=hs[a]-ad[mi[a]]+ad[mi[u]];
break;
}mi[u]=1,hs[u]=hs[a]+ad[1];
}
}g[1]=(r+1)%MOD;
for(int x=2;x<=n;x++){
ll&re=mp[hs[x]];
if(!re){
re=g[1];
for(int i=2,lm=sqrt(x);i<=lm;i++)if(x%i==0){
re+=g[i];
if(x/i^i)re+=g[x/i];
}re=re%MOD*r%MOD;
}g[x]=re;
}
for(int i=2;i<=n;i++)(g[i]+=g[i-1])%=MOD;
}
int id(ll x){return x<=B?x:B+n/x;}
int main()
{
freopen("bigben.in","r",stdin);
freopen("bigben.out","w",stdout);
n=read(),C=read(),B=sqrt(n),m=min(n,MAXN*10-50ll);
// double MUDA=clock();
if(!C)return print(1),0;
if(C==1)return print(n>2?0:1),0;
for(int i=1;i<=191;i++)ad[i]=Rand();
r=C*ksm(MOD+1-C,MOD-2,MOD)%MOD,f[id(1)]=(r+1)%MOD;
sieve(m);
for(ll d=2,x;d<=n;d=x+1){
if((x=n/(n/d))<=m){f[id(x)]=g[x];continue;}
ll&res=f[id(x=n/(n/d))]=(r+1)%MOD;
for(ll i=2,j,la;i<=x;i=la+1)
j=x/i,la=x/j,(res+=(la-i+1)%MOD*r%MOD*f[id(j)])%=MOD;
}
print(f[id(n)]*ksm(MOD+1-C,n,MOD)%MOD);
// printf("%.0f\n",clock()-MUDA);
return 0;
}