题解 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$算法$Pollard-Rho$算法
在上述推导中其实我们发现这个矩阵只与这个次数$c_i$有关,这样我们可以稍微优化下,记录一个最大的$c_i$,然后用矩阵快速幂推出这个矩阵,再将其前缀和一下,这样就可以直接查找每个$c_i$下的贡献,同时在预处理时,我们也可以把相同的$c_i$合并,计算贡献时直接一个次方上去即可。

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

发表于 2018-06-21 13:45:15 in 题解