题目背景
卷爷的刀已经抵在喉头；寒光一闪；我看到鲜血飞溅。我的身体不受控制地坠向地面……
我从现实中惊醒。“幸好在梦世界里不会发生这样的事。” 我长吁一口气。
题目描述
求
n
×
n
n\times n
n×n 矩阵
A
A
A 的行列式，其中
A
i
,
j
=
{
1
(
i
=
j
)
0
(
i
≠
j
∧
i
∣
j
)
C
(
otherwise
)
A_{i,j}=
对 998244353 998244353 998244353 取模。
数据范围与提示
n
⩽
1
0
11
n\leqslant 10^{11}
n⩽1011 。时限
6
s
\rm 6s
6s 。
简单推导 det \det det 的方式：下一行减去上一行，变成海森堡矩阵。
复杂方法：将对角线上 1 1 1 视作 C + ( 1 − C ) C+(1-C) C+(1−C)，枚举哪些位置选择了 ( 1 − C ) (1-C) (1−C)，则剩余部分是主子式。
引理：当且仅当主子式保留的行列编号 { x i } \{x_i\} {xi} 满足 x 1 ∣ x 2 ∣ x 3 ∣ ⋯ ∣ x k x_1\mid x_2\mid x_3\mid\cdots\mid x_k x1∣x2∣x3∣⋯∣xk 时，其 d e t \rm det det 是 C k C^k Ck，否则是 0 0 0 。
证明：每一行都只有一个真后缀为全零，其余是 C C C 。所以必为下三角矩阵，此时两两之间有倍数关系。
由此，我们对没选择的行列编号进行
d
p
\tt dp
dp，可以很容易地写出转移
f
(
n
)
=
C
(
1
−
C
)
n
−
1
+
∑
d
∣
n
f
(
d
)
(
1
−
C
)
n
−
d
−
1
C
f(n)=C(1-C)^{n-1}+\sum_{d\mid n}f(d)(1-C)^{n-d-1}C
f(n)=C(1−C)n−1+d∣n∑f(d)(1−C)n−d−1C
以及
a
n
s
=
(
1
−
C
)
n
+
∑
i
=
1
n
f
(
i
)
(
1
−
C
)
n
−
i
ans=(1-C)^n+\sum_{i=1}^{n}f(i)(1-C)^{n-i}
ans=(1−C)n+i=1∑nf(i)(1−C)n−i
稍微变形一下
f
(
n
)
(
1
−
C
)
−
n
=
C
1
−
C
+
∑
d
∣
n
C
1
−
C
(
f
(
d
)
(
1
−
C
)
−
d
)
f(n)(1-C)^{-n}=\frac{C}{1-C}+\sum_{d\mid n}\frac{C}{1-C}(f(d)(1-C)^{-d})
f(n)(1−C)−n=1−CC+d∣n∑1−CC(f(d)(1−C)−d)
简记
v
:
=
C
1
−
C
v:=\frac{C}{1-C}
v:=1−CC，不妨设
C
≠
1
C\ne 1
C=1 。并记
g
(
n
)
:
=
f
(
n
)
(
1
−
C
)
−
n
g(n):=f(n)(1-C)^{-n}
g(n):=f(n)(1−C)−n 。
g
(
n
)
=
v
+
∑
d
∣
n
d
≠
n
v
⋅
g
(
d
)
(1)
g(n)=v+\sum_{d\mid n}^{d\ne n}v\cdot g(d)\tag{1}
g(n)=v+d∣n∑d=nv⋅g(d)(1)
a n s = ( 1 − C ) n ∑ i = 1 n g ( i ) ans=(1-C)^n\sum_{i=1}^{n}g(i) ans=(1−C)ni=1∑ng(i)
注意到 g ( i ) g(i) g(i) 不是积性函数，我开始脑抽，打算写出关于质因子指数的 g ( i ) g(i) g(i) 的表达式。
设
n
=
∏
p
i
t
i
n=\prod p_i^{t_i}
n=∏piti，记
λ
l
:
=
∏
(
l
−
1
+
t
i
t
i
)
\lambda_l:=\prod{l-1+t_i\choose t_i}
λl:=∏(til−1+ti)
即单个质因子的指数下降方案，在长为
l
l
l 的链上。容斥以确保每一步着实有下降。
α
l
:
=
∑
i
=
0
l
−
1
(
l
−
1
i
)
(
−
1
)
i
λ
l
−
i
\alpha_l:=\sum_{i=0}^{l-1}{l-1\choose i}(-1)^i\lambda_{l-i}
αl:=i=0∑l−1(il−1)(−1)iλl−i
g
(
n
)
=
∑
i
=
1
log
n
α
i
v
i
=
∑
i
=
1
log
n
v
i
∑
j
=
0
i
−
1
(
i
−
1
j
)
(
−
1
)
j
λ
i
−
j
=
∑
i
=
1
log
n
v
i
∑
j
=
0
i
−
1
(
i
−
1
j
)
(
−
1
)
j
∏
(
i
−
j
−
1
+
t
k
t
k
)
=
∑
l
(
∑
i
=
l
+
1
log
n
v
i
(
i
−
1
i
−
l
−
1
)
(
−
1
)
i
−
l
−
1
)
∏
(
l
+
t
k
t
k
)
=
∑
l
β
l
⋅
I
l
+
1
(
n
)
I l + 1 I^{l+1} Il+1 指 ( l + 1 ) (l{+}1) (l+1) 个 I I I 的狄利克雷卷积的结果。于是可以 O ( log n ) \mathcal O(\log n) O(logn) 次筛法解决，比如杜教筛或 min25 \texttt{min25} min25 。
杜教筛不基于积性函数。而且，“因子的值求和” 不像狄利克雷卷积吗？这不像分治 FTT \textit{FTT} FTT 吗？
事实上，
(
1
)
(1)
(1) 式两边同时求前缀和立刻有
⇒
S
(
n
)
=
v
n
+
∑
i
=
2
n
v
⋅
S
(
⌊
n
/
i
⌋
)
\Rightarrow S(n)=vn+\sum_{i=2}^{n}v\cdot S(\lfloor n/i\rfloor)
⇒S(n)=vn+i=2∑nv⋅S(⌊n/i⌋)
其中 S ( n ) = ∑ i ⩽ n g ( i ) S(n)=\sum_{i\leqslant n}g(i) S(n)=∑i⩽ng(i) 。于是杜教筛吧。时间复杂度 O ( n 2 / 3 ) \mathcal O(n^{2/3}) O(n2/3) 的……吗？
注意非积性函数不能做严格意义上的线性筛。需优化。事实上，质因数的指数的可重集相同时， g g g 的值显然相同。因此我们只有很少的值需要枚举因子。但是怎么找到对应关系呢？
一个方法是 hash \texttt{hash} hash 。当然我们也可以魔改线性筛，每个数记录其最大质因子和次数，并预处理每个质数的幂。然后存 f i f_i fi 为，与 i i i 拥有相同的指数可重集的数中，最小的一个。直接求 f i f_i fi 麻烦，可以求 j j j 使得 f i = f j f_i=f_j fi=fj 。设 i i i 移除最大质因子后得到 u u u，若 f u ≠ u f_u\ne u fu=u 则 j = i u f u j=\frac{i}{u}f_u j=uifu，否则再做指数与质因数编号的判断。
最终只有 500 + 500+ 500+ 个数需暴力计算，所以复杂度总算是 O ( n 2 / 3 ) \mathcal O(n^{2/3}) O(n2/3) 了。
#include <cstdio>
#include <algorithm> // Almighty XJX yyds!!
#include <cstring> // oracle: ZXY yydBUS!!!
#include <cctype> // I hate rainybunny.
using llong = long long;
# define rep(i,a,b) for(int i=(a); i<=(b); ++i)
# define drep(i,a,b) for(int i=(a); i>=(b); --i)
# define rep0(i,a,b) for(int i=(a); i!=(b); ++i)
inline int readint(){
int a = 0, c = getchar(), f = 1;
for(; !isdigit(c); c=getchar()) if(c == '-') f = -f;
for(; isdigit(c); c=getchar()) a = a*10+(c^48);
return a*f;
}
const int LOGN = 40, MOD = 998244353;
inline llong qkpow(llong b, int q){
llong a = 1;
for(; q; q>>=1,b=b*b%MOD) if(q&1) a = a*b%MOD;
return a;
}
const int MAXN = 23544347;
bool isPrime[MAXN]; int* ppow[MAXN];
int primes[MAXN], primes_size, g[MAXN];
int fa[MAXN], big[MAXN], num[MAXN], nxt[MAXN], cnt[MAXN];
void dfs(int &res, int x, int y){
if(x == 1){ if(g+y != &res){ res = (res+g[y])%MOD; } return; }
rep(i,0,num[x]) dfs(res,nxt[x],y), y *= big[x];
}
void sieve(const int &v, int n = MAXN-1){
memset(isPrime+2, true, n-1);
fa[1] = 1, g[1] = v;
for(int i=2,&len=primes_size; i<=n; ++i){
if(isPrime[i]){
primes[++len] = big[i] = i, cnt[i] = 1;
num[i] = 1, fa[i] = 2, nxt[i] = 1;
g[i] = int(llong(v+1)*v%MOD);
for(int j=2,p=i; true; ++j,p*=i)
if(p > n/i){ ppow[i] = new int[j]; break; }
for(int j=0,p=1; p<=n/i; ++j,p*=i,ppow[i][j]=p);
}
else{
const int rt = fa[nxt[i]];
if(rt != nxt[i]) g[i] = g[fa[i] = fa[rt*(i/nxt[i])]];
else if(big[i] != primes[cnt[rt]+1]){ // not nearest
fa[i] = rt*ppow[primes[cnt[rt]+1]][num[i]];
fa[i] = fa[fa[i]], g[i] = g[fa[i]];
}
else if(rt != 1 && num[rt] < num[i]){
fa[i] = nxt[rt]*ppow[big[rt]][num[i]]
*ppow[big[i]][num[rt]]; // switch
fa[i] = fa[fa[i]], g[i] = g[fa[i]];
}
else{ // I am smallest
static int wxk = 0; ++ wxk;
fa[i] = i, dfs(g[i],i,1);
g[i] = int(llong(g[i]+1)*v%MOD);
}
}
for(int j=1; j<=len; ++j){
const int to = i*primes[j]; if(to > n) break;
isPrime[to] = false, big[to] = big[i];
if(!(i%primes[j])){
if(nxt[i] == 1) num[to] = num[i]+1, nxt[to] = 1;
else num[to] = num[i], nxt[to] = nxt[i]*primes[j];
cnt[to] = cnt[i]; break;
}
num[to] = num[i], nxt[to] = nxt[i]*primes[j], cnt[to] = cnt[i]+1;
}
}
}
int res[4650];
int getSum(const llong &n, const int &v, llong x){
if(x < MAXN) return g[x];
int &w = res[n/x]; if(~w) return w;
w = int(x%MOD); // easiest
for(llong l=2,r,val; l!=x+1; l=r){
val = x/l, r = x/val+1; // move
w = int((w+(r-l)%MOD*getSum(n,v,val))%MOD);
}
return w = int(llong(w)*v%MOD);
}
int main(){
llong n; scanf("%lld",&n); int c = readint();
if(c == 1){ puts(n <= 2 ? "1" : "0"); return 0; }
int v = int(c*qkpow(1-c+MOD,MOD-2)%MOD);
sieve(v); rep0(i,2,MAXN) g[i] = (g[i]+g[i-1])%MOD;
memset(res, -1, sizeof(res));
int ans = getSum(n,v,n)+1;
ans = int(ans*qkpow(1+MOD-c,int(n%(MOD-1)))%MOD);
printf("%d\n",ans);
return 0;
}