[acnoi2022] I have done it, but I can't

subject

Background
The rolling master's knife has reached his throat ; A flash of cold light ; I saw blood spattering . My body fell uncontrollably to the ground ……

I woke up from reality .“ Fortunately, it doesn't happen in the dream world .” I take a long breath .

Title Description
seek $n\times n$ matrix $A$ The determinant of , among
$$A_{i,j}=\{\begin{array}{ll}1& (i=j)\\ 0& (i\ne j\wedge i\mid j)\\ C& (\text{otherwise})\end{array}$$

Yes $998244353$ modulus .

Data range and tips
$n⩽10^{11}$ . time limit $6\mathrm{s}$ .

Ideas

Simple deduction $\det$ The way ： Subtract the previous line from the next line , Become Heisenberg matrix .

Complex method ： Place diagonal $1$ Regard as $C+(1-C)$, Enumerate which locations are selected $(1-C)$, Then the remaining part is the main subformula .

lemma ： If and only if the row and column numbers reserved by the primary and secondary formulas $\{x_i\}$ Satisfy $x_1\mid x_2\mid x_3\mid \cdots \mid x_k$ when , Its $\mathrm{d}\mathrm{e}\mathrm{t}$ yes $C^k$, It is $0$ .

prove ： Each line has only one true suffix of all zeros , The rest is $C$ . So it must be a lower triangular matrix , At this time, there is a multiple relationship between two pairs .

$\square$

thus , We number the unselected rows and columns $\mathtt{d}\mathtt{p}$, It is easy to write a transfer
$$f(n)=C(1-C{)}^{n-1}+\sum _{d\mid n}f(d)(1-C{)}^{n-d-1}C$$

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

Deform a little
$$f(n)(1-C{)}^{-n}=\frac{C}{1-C}+\sum _{d\mid n}\frac{C}{1-C}(f(d)(1-C{)}^{-d})$$

Brief notes $v:=\frac{C}{1-C}$, Might as well set $C\ne 1$ . Syncopation $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)$$

be aware $g(i)$ It's not an integral function , I began to think , I'm going to write about the quality factor index $g(i)$ The expression of .

Cerebral pumping

set up $n=\prod p_i^{t_i}$, remember
$${\lambda }_l:=\prod \left(\genfrac{}{}{0ex}{}{l-1+t_i}{t_i}\right)$$

That is, the exponential descent scheme of a single qualitative factor , In the long for $l$ On the chain . Inclusion and exclusion to ensure that every step really drops .
$${\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}$ finger $(l+1)$ individual $I$ The result of the Dirichlet convolution of . So you can $\mathcal{O}(\log n)$ Solve the problem by secondary screening , Such as Du Jiao sieve or $\text{min25}$ .

Positive solution

Duchenne sieves are not based on integral functions . and ,“ Sum the values of the factors ” Not like Dirichlet convolution ？ It's not like divide and conquer $\text{FTT}$ Do you ？

in fact ,$(1)$ Find the prefix and sum at the same time on both sides of the formula, and there is
$$\Rightarrow S(n)=vn+\sum _{i=2}^{n}v\cdot S(\lfloor n/i\rfloor )$$

among $S(n)={\sum }_{i⩽n}g(i)$ . So Du taught me to sift . Time complexity $\mathcal{O}(n^{2/3})$ Of …… Do you ？

Note that non integrable functions cannot be strictly linear sieves . Need to optimize . in fact , The exponents of prime factors are the same as the resettable ones ,$g$ Is obviously the same . So we only have a few values to enumerate factors . But how to find the corresponding relationship ？

One way is $\text{hash}$ . Of course, we can also change the linear sieve , Each number records its maximum prime factor and number , And preprocess the power of each prime number . Then save $f_i$ by , And $i$ Among the numbers that have the same exponentially resettable , The smallest one . Ask directly $f_i$ trouble , You can ask for $j$ bring $f_i=f_j$ . set up $i$ After removing the maximum prime factor, we get $u$, if $f_u\ne u$ be $j=\frac{i}{u}{f}_u$, Otherwise, judge the number of index and prime factor .

In the end only $500+$ The number needs to be calculated violently , So the complexity is $\mathcal{O}(n^{2/3})$ 了 .

Code

#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;
}