原文地址：http://www.cnblogs.com/GXZlegend/p/6825132.html
题目描述
Alice想要得到一个长度为n的序列，序列中的数都是不超过m的正整数，而且这n个数的和是p的倍数。Alice还希望，这n个数中，至少有一个数是质数。Alice想知道，有多少个序列满足她的要求。
输入
一行三个数，n,m,p。
1<=n<=10^9,1<=m<=2×10^7,1<=p<=100
输出
一行一个数，满足Alice的要求的序列数量，答案对20170408取模。
样例输入
3 5 3
样例输出
33
题解
矩阵乘法
至少有1个质数的方案数=总方案数-不含质数的方案数
可以先在O(m)的时间内把1~m的质数筛出来。
然后考虑：从mod p=0，到mod p=0，可以由mod p=a和mod p=p-a两个阶段组成。
可以设f[i][j]表示从mod p=i到mod p=j的方案数，不难看出这是一个矩阵，而且自乘m次就能得到答案。
所以只要处理出2种f即可。
对于每个数可以更新所有的f[i]（0≤i<n），不过这样会TLE
其实这样做没有必要，因为f都是循环出现的，只需要求f[0]即可，再据此推其余的f。
代码中，我将mod p=0当作了p来处理，其实看作0也是一样的。
#include <cstdio>
#include <cstring>
#include <algorithm>
#define mod 20170408
using namespace std;
typedef long long ll;
bool notprime[20000010];
int pri[5000010] , top;
struct matrix
{
int n , m;
ll num[110][110];
matrix()
{
n = m = 0;
memset(num , 0 , sizeof(num));
}
matrix operator*(const matrix a)const
{
matrix ans;
ans.n = n , ans.m = a.m;
int i , j , k;
for(i = 1 ; i <= ans.n ; i ++ )
for(j = 1 ; j <= ans.m ; j ++ )
for(k = 1 ; k <= m ; k ++ )
ans.num[i][j] = (ans.num[i][j] + num[i][k] * a.num[k][j]) % mod;
return ans;
}
}A , B;
matrix pow(matrix x , int y)
{
matrix ans;
int i;
ans.n = x.n , ans.m = x.m;
for(i = 1 ; i <= ans.n ; i ++ )
ans.num[i][i] = 1;
while(y)
{
if(y & 1) ans = ans * x;
x = x * x , y >>= 1;
}
return ans;
}
int main()
{
int n , m , p , i , j;
scanf("%d%d%d" , &n , &m , &p);
notprime[1] = 1;
for(i = 2 ; i <= m ; i ++ )
{
if(!notprime[i]) pri[++top] = i;
for(j = 1 ; j <= top && i * pri[j] <= m ; j ++ )
{
notprime[i * pri[j]] = 1;
if(i % pri[j] == 0) break;
}
}
A.n = A.m = B.n = B.m = p;
for(i = 1 ; i <= m ; i ++ )
{
A.num[p][(i - 1) % p + 1] ++ ;
if(notprime[i]) B.num[p][(i - 1) % p + 1] ++ ;
}
for(i = p - 1 ; i >= 1 ; i -- )
for(j = 1 ; j <= p ; j ++ )
A.num[i][j] = A.num[i + 1][j % p + 1] , B.num[i][j] = B.num[i + 1][j % p + 1];
printf("%lld
" , (pow(A , n).num[p][p] - pow(B , n).num[p][p] + mod) % mod);
return 0;
}