题解 P4714 【「数学」约数个数和】

这里给出矩阵快速幂的做法，个人认为矩乘比组合数好理解一些。 我们先考虑对于$p^m(p\text{为质数})$进行题目所描述的操作，我们可以先列出一个初始状态（第一行并不属于矩阵，只是表示第几列）： $\qquad\qquad\qquad\qquad\begin{pmatrix}0&1&2&3&\cdots&m-1&m\\0&0&0&0&\cdots&0&1\end{pmatrix}$ 表示对于每一个$p^k(k\in[0,m])$，它的贡献是多少，起始只有一个数$p^m$，它的贡献即为$1$，其它均为$0$。 然后我们进行一次求约数的操作，因为$p^m$的约数为$p^0,p^1,p^2,\dots,p^m$，所以矩阵就变成这样： $\qquad\qquad\qquad\qquad\begin{pmatrix}0&1&2&3&\cdots&m-1&m\\1&1&1&1&\cdots&1&1\end{pmatrix}$ 再进行一次求约数操作，因为对于之前产生的约数$p^k$，它产生的约数为$p^0,p^1,p^2,\dots,p^k$，所以矩阵就变成这样： $\qquad\qquad\qquad\begin{pmatrix}0&1&2&3&\cdots&m-1&m\\m&m-1&m-2&m-3&\cdots&2&1\end{pmatrix}$ 于是我们就找到了递推的公式： 设矩阵为$f[i][j]$，$i$表示进行了几次求约数操作，$j$则表示$p^j$的贡献，初始状态为$f[0][0\to m]=1$，因为第$0$次即为进行过一次操作，对于每一层$i(i>0)$，有递推式$f[i][j]=\sum\limits_{k=j}^mf[i-1][k]$，最后递推至$K$层时，显然$\sum\limits_{i=0}^mf[K][i]$即为$p^m$这个数的贡献。 而根据算术基本定理，$N=\prod\limits_{i=1}^mp_i^{c_i}$，我们就可以对于每一个$p_i^{c_i}$用上述方法计算出它的贡献，而整个数$N$的总贡献即为$\prod\limits_{i=1}^m\sum\limits_{j=0}^mf[K][j]$。 然后我们考虑用矩阵快速幂去加速这个递推，因为原矩阵为$1\times (m+1)$大小，要使得乘完后仍为$1\times (m+1)$，显然我们需要一个$(m+1)\times (m+1)$的递推矩阵： $\qquad\qquad\qquad\quad\begin{pmatrix}&0&1&2&\cdots&m-1&m\\0&1&0&0&\cdots&0&0\\1&1&1&0&\cdots&0&0\\2&1&1&1&\cdots&0&0\\\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\m-1&1&1&1&\cdots&1&0\\m&1&1&1&\cdots&1&1\end{pmatrix}$ 而因为$2^{60}=1,152,921,504,606,846,976>10^{18}$，所以$m<60$，跑的还是挺快的，于是剩下的时间就是将$N$分解质因数的时间了，如果使用朴素的$\sqrt{N}$算法显然会$T$，所以我们需要用$Miller-Rabin$算法和$Pollard-Rho$算法去分解，这里推荐个讲这两个算法的博客[$Miller-Rabin$算法](https://blog.csdn.net/sunshine_cfbsl/article/details/52425798)、[$Pollard-Rho$算法](https://blog.csdn.net/sunshine_cfbsl/article/details/52512706)。 在上述推导中其实我们发现这个矩阵只与这个次数$c_i$有关，这样我们可以稍微优化下，记录一个最大的$c_i$，然后用矩阵快速幂推出这个矩阵，再将其前缀和一下，这样就可以直接查找每个$c_i$下的贡献，同时在预处理时，我们也可以把相同的$c_i$合并，计算贡献时直接一个次方上去即可。 ```cpp #include<cstdio> #include<cstring> #include<cstdlib> #include<ctime> #include<algorithm> using namespace std; typedef long long ll; const int N=65; const int mod=998244353; int a[N],b[N],l,L,ma; //a记录不同的次方，b记录相同的次方个数,l为pr数组长度，L为b数组长度，ma为最大次方 ll pr[N<<4]; //pr为N分解出来的所以质因数 struct dd{//(m+1)*(m+1)的矩阵 ll T[N][N]; dd() { memset(T,0,sizeof(T)); } dd operator * (const dd& B)//自乘 { dd C; int i,j,k; for(i=0;i<=ma;i++) for(j=0;j<=ma;j++) for(k=0;k<=ma;k++) (C.T[i][j]+=T[i][k]*B.T[k][j]%mod)%=mod; return C; } }; struct jc{//1*(m+1)的矩阵 ll S[N]; jc() { memset(S,0,sizeof(S)); } jc operator * (const dd& B)//乘上(m+1)*(m+1)的矩阵 { jc C; int i,j; for(i=0;i<=ma;i++) for(j=0;j<=ma;j++) (C.S[i]+=S[j]*B.T[j][i]%mod)%=mod; return C; } }; jc an;//储存答案 dd tf;//递推矩阵 ll re()//快读 { ll x=0; char c=getchar(); for(;c<'0'||c>'9';c=getchar()); for(;c>='0'&&c<='9';c=getchar()) x=x*10+(c-'0'); return x; } void jc_ksm(ll y)//矩阵快速幂 { int i,j; memset(tf.T,0,sizeof(tf.T)); for(i=0;i<=ma;i++)//初始化为1 an.S[i]=1; for(i=0;i<=ma;i++)//初始化为一个形如三角形的1 for(j=0;j<=i;j++) tf.T[i][j]=1; for(;y;y>>=1,tf=tf*tf)//快速幂 if(y&1) an=an*tf; } ll n_ksm(ll x,int y)//一个普通的快速幂 { ll s=1; for(;y;y>>=1,x=x*x%mod) if(y&1) s=s*x%mod; return s; } void su()//对答案矩阵进行前缀和处理 { int i; for(i=ma;i;i--) (an.S[i-1]+=an.S[i])%=mod; } //以下为Miller和Porllard ll mu(ll x,ll y,ll M) { ll s=0; for(;y;y>>=1,x=(x+x)%M) if(y&1) s=(s+x)%M; return s; } ll ksm(ll x,ll y,ll M) { ll s=1; for(;y;y>>=1,x=mu(x,x,M)) if(y&1) s=mu(s,x,M); return s; } bool miller(ll x) { int i,j,t=0; ll y=x-1,X,Y; if(x==2) return true; for(;!(y&1);t++,y>>=1); for(j=1;j<=20;j++) { X=ksm(2+rand()%(x-2),y,x); for(i=1;i<=t;i++,X=Y) { Y=mu(X,X,x); if(Y==1&&X!=1&&X!=x-1) return false; } if(X!=1) return false; } return true; } ll gcd(ll x,ll y) { if(!y) return x; return gcd(y,x%y); } ll poll(ll x,ll c) { ll i,k=2,X=1+rand()%(x-1),Y=X,G; for(i=1;;i++) { X=(mu(X,X,x)+c)%x; G=gcd((Y-X+x)%x,x); if(G!=1&&G!=x) return G; if(X==Y) return x; if(i==k) { Y=X; k<<=1; } } } void sch(ll x,ll c) { ll k=c,p=x; if(x==1) return; if(miller(x)) { pr[++l]=x; return; } while(p>=x) p=poll(p,k--); sch(p,c); sch(x/p,c); } //以上为Miller和Porllard void ud(int k)//记录次方 { if(!b[k])//若没有出现 { a[++L]=k; b[k]=1; ma=ma>k?ma:k; } else//若出现过 b[k]++; } int main() { int i,k=0; ll n,m,s=1; n=re(); m=re(); sch(n,113);//分解质因数 sort(pr+1,pr+l+1);//排序 for(i=1;i<=l;i++)//记录 { if(i==1||pr[i]==pr[i-1]) k++; else { ud(k); k=1; } } if(k) ud(k); jc_ksm(m);//矩阵快速幂 su();//前缀和 for(i=1;i<=L;i++)//按上述进行累乘出答案 (s*=n_ksm(an.S[ma-a[i]],b[a[i]]))%=mod; printf("%lld",s); return 0; } ```