[ACNOI2022]做过也不会

题目

题目背景
卷爷的刀已经抵在喉头；寒光一闪；我看到鲜血飞溅。我的身体不受控制地坠向地面……

我从现实中惊醒。“幸好在梦世界里不会发生这样的事。” 我长吁一口气。

题目描述
求 $n\times n$ 矩阵 $A$ 的行列式，其中
$$A_{i,j}=\{\begin{array}{ll}1& (i=j)\\ 0& (i\ne j\wedge i\mid j)\\ C& (\text{otherwise})\end{array}$$

对 $998244353$ 取模。

数据范围与提示
$n⩽10^{11}$ 。时限 $6\mathrm{s}$ 。

思路

简单推导 $\det$ 的方式：下一行减去上一行，变成海森堡矩阵。

复杂方法：将对角线上 $1$ 视作 $C+(1-C)$，枚举哪些位置选择了 $(1-C)$，则剩余部分是主子式。

引理：当且仅当主子式保留的行列编号 $\{x_i\}$ 满足 $x_1\mid x_2\mid x_3\mid \cdots \mid x_k$ 时，其 $\mathrm{d}\mathrm{e}\mathrm{t}$ 是 $C^k$，否则是 $0$ 。

证明：每一行都只有一个真后缀为全零，其余是 $C$ 。所以必为下三角矩阵，此时两两之间有倍数关系。

$\square$

由此，我们对没选择的行列编号进行 $\mathtt{d}\mathtt{p}$，可以很容易地写出转移
$$f(n)=C(1-C{)}^{n-1}+\sum _{d\mid n}f(d)(1-C{)}^{n-d-1}C$$

以及
$$ans=(1-C{)}^n+\sum _{i=1}^{n}f(i)(1-C{)}^{n-i}$$

稍微变形一下
$$f(n)(1-C{)}^{-n}=\frac{C}{1-C}+\sum _{d\mid n}\frac{C}{1-C}(f(d)(1-C{)}^{-d})$$

简记 $v:=\frac{C}{1-C}$，不妨设 $C\ne 1$ 。并记 $g(n):=f(n)(1-C{)}^{-n}$ 。
$$g(n)=v+\sum _{d\mid n}^{d\ne n}v\cdot g(d)\text{\hspace{1em}(1)}$$

$$ans=(1-C{)}^n\sum _{i=1}^{n}g(i)$$

注意到 $g(i)$ 不是积性函数，我开始脑抽，打算写出关于质因子指数的 $g(i)$ 的表达式。

脑抽中

设 $n=\prod p_i^{t_i}$，记
$${\lambda }_l:=\prod \left(\genfrac{}{}{0ex}{}{l-1+t_i}{t_i}\right)$$

即单个质因子的指数下降方案，在长为 $l$ 的链上。容斥以确保每一步着实有下降。
$${\alpha }_l:=\sum _{i=0}^{l-1}\left(\genfrac{}{}{0ex}{}{l-1}{i}\right)(-1{)}^i{\lambda }_{l-i}$$

$$\begin{array}{rl}g(n)& =\sum _{i=1}^{\log n}{\alpha }_i{v}^i=\sum _{i=1}^{\log n}{v}^i\sum _{j=0}^{i-1}\left(\genfrac{}{}{0ex}{}{i-1}{j}\right)(-1{)}^j{\lambda }_{i-j}\\ & =\sum _{i=1}^{\log n}{v}^i\sum _{j=0}^{i-1}\left(\genfrac{}{}{0ex}{}{i-1}{j}\right)(-1{)}^j\prod (\genfrac{}{}{0ex}{}{i-j-1+t_k}{t_k})\\ & =\sum _l\left(\sum _{i=l+1}^{\log n}{v}^i\left(\genfrac{}{}{0ex}{}{i-1}{i-l-1}\right)(-1{)}^{i-l-1}\right)\prod \left(\genfrac{}{}{0ex}{}{l+t_k}{t_k}\right)\\ & =\sum _l{\beta }_l\cdot I^{l+1}(n)\end{array}$$

$I^{l+1}$ 指 $(l+1)$ 个 $I$ 的狄利克雷卷积的结果。于是可以 $\mathcal{O}(\log n)$ 次筛法解决，比如杜教筛或 $\text{min25}$ 。

正解

杜教筛不基于积性函数。而且，“因子的值求和” 不像狄利克雷卷积吗？这不像分治 $\text{FTT}$ 吗？

事实上，$(1)$ 式两边同时求前缀和立刻有
$$\Rightarrow S(n)=vn+\sum _{i=2}^{n}v\cdot S(\lfloor n/i\rfloor )$$

其中 $S(n)={\sum }_{i⩽n}g(i)$ 。于是杜教筛吧。时间复杂度 $\mathcal{O}(n^{2/3})$ 的……吗？

注意非积性函数不能做严格意义上的线性筛。需优化。事实上，质因数的指数的可重集相同时，$g$ 的值显然相同。因此我们只有很少的值需要枚举因子。但是怎么找到对应关系呢？

一个方法是 $\text{hash}$ 。当然我们也可以魔改线性筛，每个数记录其最大质因子和次数，并预处理每个质数的幂。然后存 $f_i$ 为，与 $i$ 拥有相同的指数可重集的数中，最小的一个。直接求 $f_i$ 麻烦，可以求 $j$ 使得 $f_i=f_j$ 。设 $i$ 移除最大质因子后得到 $u$，若 $f_u\ne u$ 则 $j=\frac{i}{u}{f}_u$，否则再做指数与质因数编号的判断。

最终只有 $500+$ 个数需暴力计算，所以复杂度总算是 $\mathcal{O}(n^{2/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;
}