bzoj4818

http://www.lydsy.com/JudgeOnline/problem.php?id=4818

矩阵快速幂+dp

首先我们来写一个dp dp[i][j]:选到第i个数，和为j，复杂度nm,不行，那么我们把j模p一下，复杂度np，还是不行。问题出在n上，那么我们要把n优化掉。

那么我们用矩阵快速幂。

首先构造列向量，dp[0],dp[p-1]

构造系数矩阵，这里我们需要总-没有质数。那么总的矩阵很好构造，当前行为i，列为j，那么我们希望(j+x)%p=i，x=(i-j)%p，那么我们构造好了，没有质数的矩阵把质数挖掉就好了。

行为i表示dp[i],也就是和%p=i，列为j表示当前选的数%p为j，那么dp[i]=tot[j]*dp[x],(j+x)%p=i,tot表示一共有多少数%p=j。

那么就构造好了。系数矩阵倍增，乘上列向量。。。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 101, M = 20000010, mod = 20170408;
struct mat {
    ll a[N][N];
} g1, g2, A1, A2;
int n, m, p;
int pri[M / 3], mark[M], x1[N], x2[N];
void Init()
{
    mark[1] = 1;
    for(int i = 2; i <= m; ++i)
    {    
        if(!mark[i]) pri[++pri[0]] = i;
        for(int j = 1; j <= pri[0] && i * pri[j] <= m; ++j)
        {
            mark[i * pri[j]] = 1;
            if(i % pri[j] == 0) break;
        }
    }
    for(int i = 1; i <= m; ++i) ++x1[i % p], x2[i % p] += mark[i];        
}
mat operator * (mat A, mat B)
{
    mat ret; memset(ret.a, 0, sizeof(ret.a));
    for(int i = 0; i < p; ++i)
         for(int j = 0; j < p; ++j)
              for(int k = 0; k < p; ++k) ret.a[i][j] = (ret.a[i][j] + A.a[i][k] * B.a[k][j]) % mod;
    return ret;
}
mat power(mat x, int t)
{
    mat ret; memset(ret.a, 0, sizeof(ret.a));
    for(int i = 0; i < p; ++i) ret.a[i][i] = 1;
    for(; t; t >>= 1, x = x * x) if(t & 1) ret = ret * x;
    return ret;
}
int main()
{
    scanf("%d%d%d", &n, &m, &p);
    Init();
    for(int i = 0; i < p; ++i)
        for(int j = 0; j < p; ++j) 
            g1.a[i][j] += x1[((i - j) % p + p) % p],
            g2.a[i][j] += x2[((i - j) % p + p) % p];
    for(int i = 0; i < p; ++i) A1.a[i][0] = x1[i], A2.a[i][0] = x2[i];
    A1 = power(g1, n - 1) * A1;
    A2 = power(g2, n - 1) * A2;
    printf("%lld\n", (A1.a[0][0] - A2.a[0][0] + mod) % mod);               
    return 0;
}