「题解」machine

machine 一道数论题目。

题目

题目描述

一句话题意:给定 \(f _ 0(n)\),有 \(m\) 组询问,每组询问给定 \(k,x\),求 \((f _ 0\ast 1^k)(x)\pmod{998244353}\)。

想看详细题意可以看原题 pdf。

数据范围

对于 \(10\%\) 的数据,\(m=0\);
对于 \(20\%\) 的数据,\(m\leq 5\);
对于另 \(30\%\) 的数据,\(f _ 0[i]\) 全都是 \(1\);
对于另 \(30\%\) 的数据,所有询问的 \(k\) 相等;
对于 \(100\%\) 的数据,\(n\leq 10^5,m\leq 1000,0\leq f _ 0[i]<998244353,k\leq 10^5,x\leq n\)。

时空限制

  • machine.cpp/in/out
  • Time limit: 4s
  • Memory limit: 512MB

题解

解法 A(人口普查)

对于 \(10\%\) 的数据,\(m=0\)。

我们直接提交一个可以通过编译的空白主函数程序即可。

解法 B(小发现)

对于 \(20\%\) 的数据,\(m\leq 5\)。

我们不妨考虑直接模拟 \(\texttt{Dirichlet}\) 卷积的过程,并且加一点点的小优化(询问按 \(k\) 排序),时间复杂度为 \(\Theta(mn\ln n\max\{k\})\),显然过不掉。

考虑 \(\texttt{Dirichlet}\) 卷积的式子。
\[f _ k(n)=\sum _ {d\mid n}f _ {k-1}(d)\]

我们不难发现,\(\texttt{Dirichlet}\) 卷积既符合交换律又符合结合律

所以我们可以通过快速幂的方式优化卷积过程,时间复杂度为 \(\Theta(mn\ln n\log _ 2\max\{k\})\)。

解法 C(小推广)

对于另 \(30\%\) 的数据,所有询问的 \(k\) 相等。

与解法 B 相同的,我们用快速幂优化卷积过程,最后的时间复杂度为 \(\Theta(n\ln n\log _ 2k+m)\)。

解法 D(臃肿的写法)

这是 std 的解法。

对于一个数 \(n\),我们不难分解它为
\[n=\prod _ {p _ i\in\mathbb{P} _ n}p _ i^{c _ i}\]

对于一个数 \(n\),我们根据分解的结果定义一个数列
\[c _ 1,c _ 2,\cdots,c _ t\]

其中有
\[c _ 1\leq c _ 2\leq c _ 3\leq\cdots\leq c _ t\]

我们规定,如果两个数对应的数列相同,我们就称这两个数是乘法同构的

不难观察到(打表即可),对于小于等于 \(10^5\)(本题 \(n\) 的范围)的数中,根据乘法同构分类,本质不同的数仅有 \(160\) 个(题解里面写的是 \(165\) 个,显然是不认真思考,抄别人代码写的题解)。

于是我们考虑用乘法同构矩阵优化 \(\texttt{Dirichlet}\) 卷积过程

最后的时间复杂度为
\[\Theta(n\log^2 _ 2n+n\ln n+|S|^3\log _ 2k+m(|S|^2\log _ 2k+\sqrt{x}))\]

奇慢无比。

代码

std 的做法太慢了代码还长

#include<bits/stdc++.h>
using namespace std;
#define reg register
typedef long long ll;
#define MOD 998244353
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++)
static char buf[100000],*p1=buf,*p2=buf;
inline int read(void){
	reg char ch=getchar();
	reg int res=0;
	while(ch<'0'||'9'<ch)ch=getchar();
	while('0'<=ch&&ch<='9')res=10*res+ch-'0',ch=getchar();
	return res;
}

inline int add(reg int a,reg int b){
	reg int sum=a+b;
	return sum>=MOD?sum-MOD:sum;
}

const int MAXN=100000+5;
const int MAXLOG2K=17+1;

int n,m;
int f[MAXN];
int col[MAXN];
int link[MAXN];
map<vector<int>,int> S;

bool vis[MAXN];
int tot,prime[MAXN];
int from[MAXN];

inline void Init(reg int n){
	vis[1]=true;
	from[1]=1;
	for(reg int i=2;i<=n;++i){
		if(!vis[i]){
			prime[++tot]=i;
			from[i]=i;
		}
		for(reg int j=1;j<=tot&&i*prime[j]<=n;++j){
			vis[i*prime[j]]=true;
			from[i*prime[j]]=prime[j];
			if(i%prime[j]==0)
				break;
		}
	}
	return;
}

const int MAXCNT=160+1+5;

struct Matrix{
	int n,m;
	int unit[MAXCNT][MAXCNT];
	inline Matrix(void){
		n=m=0;
		memset(unit,0,sizeof(unit));
		return;
	}
	inline int* operator[](reg int i){
		return unit[i];
	}
	inline Matrix operator*(const Matrix &a)const{
		Matrix res;
		res.n=n,res.m=a.m;
		reg int r;
		for(reg int i=1;i<=n;++i)
			for(reg int k=1;k<=m;++k){
				r=unit[i][k];
				for(reg int j=1;j<=a.m;++j)
					res[i][j]=add(res[i][j],1ll*r*a.unit[k][j]%MOD);
			}
		return res;
	}
};

Matrix base[MAXLOG2K];

int main(void){
	//freopen("machine.in","r",stdin);
	//freopen("machine.out","w",stdout);
	n=read();
	for(reg int i=1;i<=n;++i)
		f[i]=read();
	Init(1e5);
	reg int cnt=0;
	col[++cnt]=1;
	link[1]=1;
	for(reg int i=2;i<=n;++i){
		vector<int> V;
		reg int tmp=i;
		while(tmp!=1){
			reg int p=from[tmp];
			int cnt=0;
			while(tmp%p==0)
				tmp/=p,++cnt;
			V.push_back(cnt);
		}
		sort(V.begin(),V.end());
		if(S.count(V))
			col[i]=S[V];
		else{
			col[i]=S[V]=++cnt;
			link[i]=1;
		}
	}
	base[0].n=base[0].m=cnt;
	for(reg int i=1;i<=n;++i)
		for(reg int k=1;i*k<=n;++k)
			base[0][col[i]][col[i*k]]+=link[i*k];
	for(reg int i=1;i<MAXLOG2K;++i)
		base[i]=base[i-1]*base[i-1];
	m=read();
	while(m--){
		static int k,x;
		k=read(),x=read();
		Matrix A;
		A.n=1,A.m=cnt;
		A[1][1]=1;
		for(reg int i=0;i<MAXLOG2K;++i)
			if((k>>i)&1)
				A=A*base[i];
		reg int ans=0;
		for(reg int i=1;i*i<=x;++i)
			if(x%i==0){
				ans=add(ans,1ll*A[1][col[x/i]]*f[i]%MOD);
				if(i*i!=x)
					ans=add(ans,1ll*A[1][col[i]]*f[x/i]%MOD);
			}
		printf("%d\n",ans);
	}
	return 0;
}

解法 E(人类智慧)

考虑递推式
\[
\begin{aligned}
f _ k(n)&=\sum _ {d _ 1\mid n}f _ {k-1}(d _ 1)\\
&=\sum _ {d _ 1\mid n}\sum _ {d _ 2\mid d _ 1}f _ {k-2}(d _ 2)\\
&=\sum _ {d _ 1\mid n}\sum _ {d _ 2\mid d _ 1}\sum _ {d _ 3\mid d _ 2}f _ {k-3}(d _ 3)\\
&=\sum _ {d _ 1\mid n}\sum _ {d _ 2\mid d _ 1} \cdots \sum _ {d _ {i+1}\mid d _ i}\cdots\sum _ {d _ k\mid d _ {k-1}}f _ 0(d _ k)
\end{aligned}
\]

我们不妨变换一下形式
\[
\begin{aligned}
f _ k(n)&=\sum _ {d _ 1\mid n}\sum _ {d _ 2\mid d _ 1} \cdots \sum _ {d _ {i+1}\mid d _ i}\cdots\sum _ {d _ k\mid d _ {k-1}}f _ 0(d _ k)\\
&=\sum _ {d _ k\mid d _ {k-1}}f _ 0(d _ k)[d _ k\mid d _ {k-1}\mid\cdots\mid d _ 2\mid d _ 1\mid n]\\
&=\sum _ {d _ k\mid n}f _ 0(d _ k)[d _ k\mid d _ {k-1}\mid\cdots\mid d _ 2\mid d _ 1\mid n]
\end{aligned}
\]

此时我们注意到后面的特征函数
\[
[d _ k\mid d _ {k-1}\mid\cdots\mid d _ 2\mid d _ 1\mid n]
\]

我们把上面的 \(d\) 看作分解后得到的集合,则
\[
\mathbb{P} _ {d _ k}\subseteq\mathbb{P} _ {d _ {k-1}}\subseteq\mathbb{P} _ {d _ {k-2}}\subseteq\cdots\subseteq\mathbb{P} _ {d _ 2}\subseteq\mathbb{P} _ {d _ 1}\subseteq\mathbb{P} _ n
\]

不难发现,我们本质上就是要求满足上述子集链关系的方案数,方案数与 \(f _ 0(d _ k)\) 的乘积即为对答案的贡献。

我们考虑把所有集合都消去一个 \(\mathbb{P} _ {d _ k}\),这对答案没有影响,然后得到
\[
\begin{matrix}
\mathbb{P} _ {d _ k}&\subseteq&\mathbb{P} _ {d _ {k-1}}&\subseteq&\cdots&\subseteq&\mathbb{P} _ {d _ 1}&\subseteq&\mathbb{P} _ n\\
\downarrow&&\downarrow&&\downarrow&&\downarrow&&\downarrow\\
\varnothing&\subseteq&\mathbb{S} _ 1&\subseteq&\cdots&\subseteq&\mathbb{S} _ {k-1}&\subseteq&(\mathbb{P} _ n-\mathbb{P} _ {d _ k})
\end{matrix}
\]

\[
\varnothing\subseteq\mathbb{S} _ 1\subseteq\mathbb{S} _ 2\subseteq\cdots\subseteq\mathbb{S} _ {k-2}\subseteq\mathbb{S} _ {k-1}\subseteq(\mathbb{P} _ n-\mathbb{P} _ {d _ k})
\]

考虑子集符号象征着集合之间的差(变化),我们可以分开质因子讨论。

上述链中共有 \(k\) 个子集符号。

对于一个质因子 \(p _ i\),它在 \((\mathbb{P} _ n-\mathbb{P} _ {d _ k})\) 中的出现次数为 \(c _ i\),那么它单独的方案数目可以转化为:

\(c _ i\) 个相同小球放入 \(k\) 个不同盒子中,允许有的盒子空着的方案数。

用插板法可以求得为 \(\binom{k+c _ i-1}{k-1}=\binom{k+c _ i-1}{c _ i}\)。

再用乘法原理对质因子对方案数的贡献进行合并,直接乘起来即可。

最后得出答案的表达式
\[
f _ k(x)=\sum _ {d\mid x}f _ 0(d)\prod _ {c _ i\in(\mathbb{P} _ x-\mathbb{P} _ d)}\binom{k+c _ i-1}{c _ i}\pmod{998244353}
\]

这种做法快到飞起,时间复杂度为 \(\Theta(n+m\sum _ {d\mid x}\sqrt{d})\),其中 \(\Theta(n)\) 的复杂度用于预处理阶乘及其逆元,实际上的复杂度跑不满,实际测试趋近于 \(\Theta(n+m\log _ 2m)\)。

代码

时间效率是解法 D 的 200 倍(平均一个点 0.01s),并且代码还短了一半。

#include<bits/stdc++.h>
using namespace std;
#define reg register
typedef long long ll;
#define MOD 998244353
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++)
static char buf[100000],*p1=buf,*p2=buf;
inline int read(void){
	reg char ch=getchar();
	reg int res=0;
	while(ch<'0'||'9'<ch)ch=getchar();
	while('0'<=ch&&ch<='9')res=10*res+ch-'0',ch=getchar();
	return res;
}

const int MAXN=100000+5;
const int MAXS=500000+5;

inline int add(reg int a,reg int b){
	reg int sum=a+b;
	return sum>=MOD?sum-MOD:sum;
}

int n;
int f[MAXN];
int fac[MAXS],invfac[MAXS];

inline int pow(reg int x,reg int exp){
	reg int res=1;
	while(exp){
		if(exp&1)
			res=1ll*res*x%MOD;
		x=1ll*x*x%MOD;
		exp>>=1;
	}
	return res;
}

inline void Init(reg int n){
	fac[1]=invfac[1]=1;
	for(reg int i=2;i<=n;++i)
		fac[i]=1ll*fac[i-1]*i%MOD;
	invfac[n]=pow(fac[n],MOD-2);
	for(reg int i=n-1;i>=2;--i)
		invfac[i]=1ll*invfac[i+1]*(i+1ll)%MOD;
	return;
}

inline int C(reg int n,reg int m){
	return 1ll*fac[n]*invfac[m]%MOD*invfac[n-m]%MOD;
}

inline int Solve(reg int x,reg int k){
	reg int res=1;
	for(reg int i=2;i*i<=x;++i)
		if(x%i==0){
			reg int cnt=0;
			while(x%i==0)
				++cnt,x/=i;
			res=1ll*res*C(k+cnt-1,cnt)%MOD;
		}
	if(x>1)
		res=1ll*res*k%MOD;
	return res;
}

int main(void){
	Init(5e5);
	n=read();
	for(reg int i=1;i<=n;++i)
		f[i]=read();
	reg int m=read();
	while(m--){
		static int k,x;
		k=read(),x=read();
		reg int ans=0;
		for(reg int i=1;i*i<=x;++i)
			if(x%i==0){
				ans=add(ans,1ll*f[i]*Solve(x/i,k)%MOD);
				if(i*i!=x)
					ans=add(ans,1ll*f[x/i]*Solve(i,k)%MOD);
			}
		printf("%d\n",ans);
	}
	return 0;
}