题解 P5616 【[MtOI2019]恶魔之树】

· · 题解

可能更好的体验

发现 s_i 最大只有 300，而 \sqrt {300} 以内的质数只有 2,3,5,7,11,13,17，也就是说，剩下的质因子在每个数中至多出现一个，且次数最多为 1。

然后，2,3,5,7,11,13,17 的最大次数分别为 8,5,3,2,2,2,2，组合一下，总方案数只有 17496 个，是个挺小的状态数。

所以我们先不考虑那些大于 \sqrt{300} 的质数产生的贡献。

但是 n 特别大，如果进行 n 次转移，那仍然接受不了。

注意只有 300 种不同的数，如果一个数出现 c 次，那么这个数对最小公倍数产生贡献的方案数就为 2^{c}-1。

那么我们对 300 以内每个最大质因子不超过 17 的数进行考虑。

以下数组都忽略第一维的滚动维。

令 dp[a1][a2][a3][a4][a5][a6][a7] 表示最小公倍数为 2^{a1}3^{a2}5^{a3}7^{a4}11^{a5}13^{a6}17^{a7} 的方案数。

转移时枚举a1,a2,a3,a4,a5,a6,a7，计算出新的最小公倍数 b1,b2,b3,b4,b5,b6,b7，然后直接转移即可。

接下来考虑大于 17 的质因数。我们把最大质因子相同的数放在一起考虑。然后考虑下一个质因子的时候，原来的质因子的有无就不会对最小公倍数产生影响。

这里不能直接计算方案数，因为不知道最大质因子的出现情况。所以直接统计总贡献。

令 F[b][a1][a2][a3][a4][a5][a6][a7] 表示最小公倍数的前 7 种质因数的情况为 2^{a1}3^{a2}5^{a3}7^{a4}11^{a5}13^{a6}17^{a7} 时，当前考虑的大质因子次数为 b 的总和。

转移时，我们计算出新的最小公倍数与原来的最小公倍数差的那些因子的乘积 \Delta，表示原来的那些所有方案都要乘上 \Delta 转移到新的状态。

这个和上面的转移也是类似的。

我们处理完一个最大质因子后，需要把 b=1 的所有情况都加给 b=0 的情况上，表示最大质因子换了一个。

那么总共的运算次数是 O(17496\times 300)。在 1 秒内跑完绰绰有余。

Code：

#include<iostream>
#include<vector>
#include<cstring>
#include<algorithm>
#define rep(i,a,b)for(int i=(a);i<=(b);++i)
using namespace std;
typedef long long LL;
//2:8 3:5 5:3 7:2 11:2 13:2 17:2
//4 3 2 2 2 2 2
int dp[2][9][6][4][3][3][3][3],F[2][2][9][6][4][3][3][3][3];
int n,md;
int cnt[305];
vector<pair<int,int> >vc[305];
vector<int>pc[305];
inline int pow(int a,int b){
    int ret=1;
    for(;b;b>>=1,a=(LL)a*a%md)if(b&1)ret=(LL)ret*a%md;
    return ret;
}
void solve(int n,vector<pair<int,int> >&vec){
    for(int i=2;i*i<=n;++i)if(n%i==0){
        int cnt=0;
        while(n%i==0)++cnt,n/=i;
        vec.emplace_back(i,cnt);
    }
    if(n>1)vec.emplace_back(n,1);
}
int main(){
    ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
    cin>>n>>md;
    for(int i=1;i<=n;++i){
        int x;
        cin>>x;
        ++cnt[x];
    }
    ********dp=1;
    int cur=0;
    for(int i=1;i<=300;++i)if(cnt[i]){
        solve(i,vc[i]);
        if(i==1||vc[i].back().first<=17){
            int fs=pow(2,cnt[i])-1;
            cur^=1;
            memcpy(dp[cur],dp[cur^1],sizeof*dp);
            rep(a1,0,8)rep(a2,0,5)rep(a3,0,3)rep(a4,0,2)rep(a5,0,2)
            rep(a6,0,2)rep(a7,0,2){
                int s=dp[cur^1][a1][a2][a3][a4][a5][a6][a7];
                if(!s)continue;
                int b1=0,b2=0,b3=0,b4=0,b5=0,b6=0,b7=0;
                for(auto j:vc[i]){
                    switch(j.first){
                        case 2:b1=j.second;break;
                        case 3:b2=j.second;break;
                        case 5:b3=j.second;break;
                        case 7:b4=j.second;break;
                        case 11:b5=j.second;break;
                        case 13:b6=j.second;break;
                        case 17:b7=j.second;break;
                    }
                }
                b1=max(b1,a1);
                b2=max(b2,a2);
                b3=max(b3,a3);
                b4=max(b4,a4);
                b5=max(b5,a5);
                b6=max(b6,a6);
                b7=max(b7,a7);
                int t=pow(2,b1)*pow(3,b2)%md*pow(5,b3)%md*pow(7,b4)%md*pow(11,b5)%md*pow(13,b6)%md
                *pow(17,b7)%md;
                int&to=dp[cur][b1][b2][b3][b4][b5][b6][b7];
                to=(to+(LL)s*fs)%md;
            }
        }else
        pc[vc[i].back().first].push_back(i);
    }
    int ans=0;
    rep(a1,0,8)rep(a2,0,5)rep(a3,0,3)rep(a4,0,2)rep(a5,0,2)
    rep(a6,0,2)rep(a7,0,2)
    F[0][0][a1][a2][a3][a4][a5][a6][a7]=dp[cur][a1][a2][a3][a4][a5][a6][a7]*((LL)
    pow(2,a1)*pow(3,a2)%md*pow(5,a3)%md*pow(7,a4)%md*pow(11,a5)%md*pow(13,a6)%md
                *pow(17,a7)%md)%md;
    cur=0;
    for(int p=2;p<=300;++p){
        for(int i:pc[p]){
            int fs=pow(2,cnt[i])-1;
            cur^=1;
            memcpy(F[cur],F[cur^1],sizeof*F);
            rep(ff,0,1)rep(a1,0,8)rep(a2,0,5)rep(a3,0,3)rep(a4,0,2)rep(a5,0,2)
            rep(a6,0,2)rep(a7,0,2){
                int s=F[cur^1][ff][a1][a2][a3][a4][a5][a6][a7];
                if(!s)continue;
                int b1=0,b2=0,b3=0,b4=0,b5=0,b6=0,b7=0;
                for(auto j:vc[i]){
                    switch(j.first){
                        case 2:b1=j.second;break;
                        case 3:b2=j.second;break;
                        case 5:b3=j.second;break;
                        case 7:b4=j.second;break;
                        case 11:b5=j.second;break;
                        case 13:b6=j.second;break;
                        case 17:b7=j.second;break;
                    }
                }
                int dlt=
                pow(2,max(b1-a1,0))*
                pow(3,max(b2-a2,0))%md*
                pow(5,max(b3-a3,0))%md*
                pow(7,max(b4-a4,0))%md*
                pow(11,max(b5-a5,0))%md*
                pow(13,max(b6-a6,0))%md*
                pow(17,max(b7-a7,0))%md*pow(p,ff^1)%md;
                b1=max(b1,a1);
                b2=max(b2,a2);
                b3=max(b3,a3);
                b4=max(b4,a4);
                b5=max(b5,a5);
                b6=max(b6,a6);
                b7=max(b7,a7);
                int&to=F[cur][1][b1][b2][b3][b4][b5][b6][b7];
                to=(to+s*(LL)dlt%md*fs)%md;
            }
        }
        rep(a1,0,8)rep(a2,0,5)rep(a3,0,3)rep(a4,0,2)rep(a5,0,2)
        rep(a6,0,2)rep(a7,0,2){
            int&f=F[cur][0][a1][a2][a3][a4][a5][a6][a7];
            int&g=F[cur][1][a1][a2][a3][a4][a5][a6][a7];
            f=(f+g)%md;
            g=0;
        }
    }
    rep(a1,0,8)rep(a2,0,5)rep(a3,0,3)rep(a4,0,2)rep(a5,0,2)
    rep(a6,0,2)rep(a7,0,2)
    ans=(ans+F[cur][0][a1][a2][a3][a4][a5][a6][a7])%md;
    cout<<ans<<endl;
    return 0;
}