标签：oid 复杂 lin display 需要 type pow line ret
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; }
标签：oid 复杂 lin display 需要 type pow line ret
原文地址：http://www.cnblogs.com/19992147orz/p/6828496.html
