ZYS的黄金——解题报告

题目链接：http://www.ycoj.cf/problem/14
题目大意：n个大臣进贡的黄金都是[1，m]内的正整数。若后一个大臣的黄金数是后一个人的倍数，那么这个人就会受罚。因此他希望你能帮他求出有多少种进贡黄金的方案，不存在一个大臣进贡的黄金恰好是前一个大臣进贡黄金的倍数。两种方案不同当且仅当存在x使得第x个进贡的大臣在两种方案中进贡了的黄金数量不同。
题目分析：
一：如果用f[i][j]表示，第i个大臣进贡j的黄金的方案数。那么f[i][j]=∑f[i-1][k]（j%k!=0），其实可以转化为：f[i][j]=∑f[i-1][k]（m>=k>=1）-∑f[i-1][k]（j%k==0）。
二：由于所有的质数在n==1时，f的值为1，那么因为f[i][j]只与j的因数有关，那么我们将一个数分解质因数后可以发现，如果质因数的指数从大到小排序后，两个数的质因数指数分别相等，那么他们的f值是相同的，因为都是由同样数量质数转移而来的。所以我们只需要处理好质数的组合情况就行了。
三：分析时间复杂度后我们发现，只有当时间复杂度在O(nlogm)左右时才可以通过题目。我们通过打表可以发现，指数组合仅只有不到200种可能性。那么就可以通过矩阵二分快速幂来实现这个过程，降低时间复杂度，最后时间复杂度大致为O(200^3logn)，可以通过题目。
四：在处理指数组合的时候可以用哈希表来处理。
解题过程：
一：用哈希表处理处每个指数组合，并记录每个组合出现的次数。
二：根据1~m的数之间的整除关系来构造矩阵，矩阵i,j表示组合i–>组合j的转移方案数。
三：利用矩阵二分快速幂迅速得到答案。
最后附上AC代码：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;
typedef int ll;
ll N,m,P;
ll d[300010],pos[300010];
ll type[210],tot=0;
ll Pri[]={2,3,5,7,11,13,17,19},tp[300010];
void getMax()
{
    ll count1=1;
    for(ll i=2;i<=m;++i)
    {
        if(!d[i])
            d[i]=tp[++count1]=i;
        ll u=min(m/i,d[i]);
        for(ll j=1;j<=count1 && tp[j]<=u;++j)
            d[i*tp[j]]=tp[j];
    }
}
bool cmp(ll a,ll b)
{
    return a>b;
}
ll mapping(ll x)
{
    ll C[20],T=0;
    while(x!=1)
    {
        ll p=d[x];
        C[T]=0;
        while(d[x]==p)
            x/=p,++C[T];
        ++T;
    }
    sort(C,C+T,cmp);
    ll ans=1,p,c;
    while(T--)
        for(p=Pri[T],c=C[T];c--;)
            ans*=p;
    if(pos[ans]==-1)
    {
        pos[ans]=tot++;
        type[pos[ans]]=ans;
    }
    return pos[ans];
}
ll cnt[210];
struct matrix
{
    ll n;
    ll a[210][210];
};
matrix matrix_mul(matrix A, matrix B)
{
    matrix ret;
    ret.n = A.n;
    for(ll i=0;i<ret.n;i++)
        for(ll j=0;j<ret.n;j++)
            ret.a[i][j]=0;
    for(ll i=0;i<ret.n;i++)
        for(ll j=0;j<ret.n;j++)
            for(ll k=0;k<A.n;k++)
                ret.a[i][j]=(ret.a[i][j]+A.a[i][k]*(long long)B.a[k][j]%P)%P;
    return ret;
}
matrix unit(ll n)
{
    matrix ret;
    ret.n=n;
    for(ll i=0;i<n;i++)
    {
        for(ll j=0;j<n;j++)
        {
            if(i==j)
                ret.a[i][j]=1;
            else
                ret.a[i][j]=0;
        }
    }
    return ret;
}
void matrix_pow(matrix A,ll p,matrix &V)
{
    while(p)
    {
        if(p&1)
            V=matrix_mul(V,A);
        A=matrix_mul(A,A);
        p>>=1;
    }
}
ll num[100010];
int main()
{
    memset(pos,-1,sizeof(pos));
    scanf("%d%d%d",&N,&m,&P);
    getMax();
    for(ll i=1;i<=m;i++)
    {
        num[i]=mapping(i);
        cnt[num[i]]++;
    }
    matrix M,V;
    M.n=tot;
    for(ll i=0;i<tot;++i)
        for(ll j=0;j<tot;++j)
            M.a[i][j]=cnt[i];
    V.n=tot;
    for(ll i=0;i<tot;i++)
        V.a[0][i]=1;
    for(ll i=0;i<tot;i++)
    {
        ll x=type[i];
        for(ll f=1;f*f<=x;++f)
        {
            if(x%f==0)
            {
                if(num[f]!=0)
                    --M.a[num[f]][i];
                else
                {
                    num[f]=mapping(f);
                    --M.a[num[f]][i];
                }
                if(f*f!=x)
                {
                    if(num[x/f]!=0)
                        --M.a[num[x/f]][i];
                    else
                    {
                        num[x/f]=mapping(f);
                        --M.a[num[x/f]][i];
                    }
                }
            }
        }
    }
    N--;
    matrix_pow(M,N,V);
    ll A=0;
    for(ll i=0;i<tot;i++)
        A=(A+cnt[i]*(long long)V.a[0][i])%P;
    printf("%d\n",A);

    return 0;
}