LOj #2002. 「SDOI2017」序列计数（容斥+dp+矩阵快速幂）

Dear 丶 2022-06-10 12:45 183阅读 0赞

题目链接：  
[LOj 2002][]

题意：  
要求得到一个长度为 n 的序列，序列中的数都是不超过 m 的正整数，而且这 n个数的和是 p 的倍数。这n个数中，至少有一个数是质数。问你有多少个序列满足要求。

题解：  
根据容斥原理，因为题目要求至少有一个数是素数，用所有方案减去不含质数的方案就是答案。  
设 dp\[i\]\[j\]表示序列前 i 个数模 p 的余数为 j 时的方案数。dp\[i\]\[j\]=∑dp\[i−1\]\[(j−k)modp\]。  
如果不考虑包含质数的情况，因为p很小，我们不妨构建p∗p的矩阵，然后对于(i,j)填有多少个数x使得(i\+x)%p=j，最后矩阵快速幂一下就行了。对于不含质数的情况，我们也同样构建这样一个矩阵，然后矩阵快速幂一下就可以了。

AC代码：

#include<bits/stdc++.h>
    using namespace std;
    typedef long long ll;
    const int mod = 20170408;
    const int sz = 110;
    int n,m,p;
    struct matrix
    {
        int ma[sz][sz];
        matrix friend operator *(const matrix a,const matrix b)
        {
            matrix ret;
            memset(ret.ma,0,sizeof(ret.ma));
            for(int i=0;i<p;i++){
                for(int j=0;j<p;j++)
                {
                    for(int k=0;k<p;k++)
                    {
                        ret.ma[i][j] = (ret.ma[i][j] + 1LL * a.ma[i][k] * b.ma[k][j]%mod)%mod;
                    }
                }
            }
            return ret;
        }
    }A,ans;
    
    
    matrix multipow(matrix x, int k)
    {
        memset(ans.ma,0,sizeof(ans.ma));
        for(int i=0;i<p;i++){
            ans.ma[i][i] = 1;
        }
        while(k)
        {
            if(k&1)ans=ans*x;
            k>>=1;
            x=x*x;
        } 
        return ans;
    }
    
    bool isprime[20000010];
    int prime[2000010];
    
    void sieve()
    {
        isprime[0] = true; isprime[1] = true;
        for(int i=2,cnt=0;i<=m;i++){
            if(isprime[i]==false) prime[++cnt] = i;
            for(int j=1;j<=cnt&&i*prime[j]<=m;j++)
            {
                isprime[i*prime[j]] = true;
                if(i%prime[j]==0)break;
            }
        }
    }
    
    int dp[110];
    
    int solve1()
    {
        for (int i=1;i<=m;i++)  
        {
            dp[i%p]++;
        }
        for (int j=1;j<=m;j++)
        {   
            A.ma[(-j%p+p)%p][0]++;
        }
        for (int i=1;i<p;i++)
        {
            for (int j=0;j<p;j++)
            {
                A.ma[j][i] = A.ma[(j-1+p)%p][i-1];
            }         
        }
    //  for(int i=1;i<=m;i++)
    //  {
        //  for(int j=1;j<=m;j++){
        //      printf("%d ",A.ma[i][j]);
        //  }
        //  cout<<endl;
    //  }
        A=multipow(A,n-1);
        int ans=0;
        for (int i=0;i<p;i++) 
        {
            ans = (ans + 1LL * dp[i] * A.ma[i][0] % mod ) % mod;
        }
        //cout<<"ans1="<<ans<<endl;
        return ans;
    }
    
    int solve2()
    {
        memset(dp,0,sizeof(dp));
        for (int i=1;i<=m;i++)
        {
            if (isprime[i]) 
            {
                dp[i%p]++;
            }
        }
    
        memset(A.ma,0,sizeof(A.ma));
    
        for (int j=1;j<=m;j++)
        {
            if (isprime[j]) 
            {
                A.ma[(-j%p+p)%p][0]++;
            }
        }
    
        for (int i=1;i<p;i++)
        {
            for (int j=0;j<p;j++)
            {
                A.ma[j][i] = A.ma[(j-1+p)%p][i-1];
            }           
        }
        //for(int i=1;i<=m;i++)
    //  {
        //  for(int j=1;j<=m;j++){
        //      printf("%d ",A.ma[i][j]);
        //  }
        //  cout<<endl;
    //  }       
        A=multipow(A,n-1);
        int ans=0;
        for (int i=0;i<p;i++) 
        {
            ans = (ans + 1LL * dp[i] * A.ma[i][0] % mod ) % mod;
        }
       // cout<<"ans2="<<ans<<endl;
        return ans;
    }
    /*
    59 74 92
    12897479
    */
    int main()
    {
        //freopen("in.txt","r",stdin);
        scanf("%d%d%d",&n,&m,&p);
        sieve();
        int x=solve1();
        int y=solve2();
        printf("%d\n",(x-y+mod)%mod);
        return 0; 
     }

[LOj 2002]: https://loj.ac/problem/2002