[2022 national tournament simulation] BigBen -- determinant, Du Jiao sieve

Good question ！ Unfortunately, there is no link ！

Title Description

Answer key

The calculation method is the same as the positive solution , But the Heisenberg matrix is not needed in the derivation process .

Consider replacing all diagonal elements with polynomials $C+x$, Then the determinant of this matrix becomes a polynomial , You can substitute $x=1-C$ Find out the answer .
Combine the selected $x$ You can get ：
$$\det (A)=\sum _{S\subset \{1,2,...,n\}}\det (B_S)x^{n-|S|}$$ among $B_S$ It's a $A$ Replace the diagonal of with $C$ Then limit it to $S$ The matrix on this subset . such as ：

Observe $B_S$ You can find , If $S$ The elements in can be divided by two （ That is, after sorting , Divide the previous by the next ）, that $B_S$ It's a lower triangular matrix , The determinant is $C^{|S|}$;
Otherwise, if you can't divide by two , Then there must be $\det (B_S)=0$. This can be proved recursively ：
If $S$ There are two factors in which the number is not any other number , Then the lines of both numbers are all $C$, Then the matrix must not have the rank ;
conversely ,$S$ The largest number in must be a multiple of all numbers , Then the last row of the column is $C$, So you can delete the last row and column of the submatrix recursively . It is easy to find that the above rank is not satisfied in the end .

We make $r=\frac{C}{1-C}$, Then the problem is to find a set that satisfies the division of two elements $S$ And ：
$$Ans=(1-C{)}^n\sum _{Selementplaintwotwowholeexcept}{r}^{|S|}$$
Make $g_x={\sum }_{\max \{i|i\in S\}=x}{r}^{|S|}$（ Specially ,$g_1$ Express $S=\{1\}$ And empty sets ）, Then you can enumerate the set without the maximum value , Get the transfer formula
$$g_i=r\sum _{j|i}{g}_j$$ And initial value $g_1=r+1$.
Our aim is to find $g_x$ And , And this is just a formula that looks like a Dirichlet convolution and cannot be used to Min_25 sieve , Need to find another way .

Make $s(n)={\sum }_{i=1}^n{g}_i$, We can actually enumerate directly $S$ The state transition is obtained by the multiple relationship between the maximum value and the second maximum value in ：
$$s(n)=1+r+\sum _{i=2}^{n}s(\lfloor \frac{n}{i}\rfloor )$$ So we get a direct radical divide and conquer enumeration .

However, the complexity of this approach is $O(n^{\frac{3}{4}})$, Can only pass 80 branch , Also optimize . Here we can follow the thought of Du Jiao sieve , Preprocess first $g_i$ Before $n^{\frac{2}{3}}$ Find the prefix and the value , Then the rest of the radical divide and conquer enumeration , In this way, the complexity of radical divide and conquer can be reduced to $O(n^{\frac{2}{3}})$.

Now the problem is , Yes? $O(n)$ Find out linearly $g_i$. consider $x$ The only decomposition of , We find that if the number of times of all prime factors is the same , that $g_x$ equal （ Because this is essentially a path count going to multiples , So all prime factors are equivalent ）. We find the hash value of all elements by linear sieve , Then, if the hash values are equal, we only need to find them once .

The magic is , There are very few kinds of repeatable sets , Less than $\sqrt{LinesexsieveFanaround}$, Therefore, we select the smallest number in the same hash value and directly enumerate the factors to find .

The final total complexity is $O(n^{\frac{2}{3}})$.

Code

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