题目:
http://poj.openjudge.cn/practice/1058/
http://www.lydsy.com/JudgeOnline/problem.php?id=3328
题意:
两题的题面不太相同,抽象后题意如下。
给出一个
M×M
的矩阵
G
,求
其中 P 是质数,并且
Guideposts: 1≤N≤1018,1≤M≤5,1≤K≤5000,1≤P≤109
PYXFIB: 1≤N≤1018,M=2,1≤K≤2⋅104,1≤P≤109
题解:
不考虑
[K|i]
的话,有
∑Ni=0(Ni)Gis,t=∑Ni=0(Ni)(IN−iGi)s,t=(I+G)Ns,t
。
而
[K|i]
可以利用单位复根进行构造,显然有
∑0≤j<KωijKK=[K|i]
则原式可化为
∑i=0N∑j=0K−1ωijKK(Ni)Gis,t=1K∑j=0K−1∑i=0N(Ni)(ωjKG)is,t=1K∑j=0K−1∑i=0N(Ni)(IN−i(ωjKG)i)s,t=1K∑j=0K−1(I+ωjKG)Ns,t
因此可以得到一个 O(M3KlogN) 的做法,但是直接使用单位复根,会导致精度过低。
由于 P 是质数,并且
代码:
#include <cstdio>
#include <cstring>
typedef long long LL;
const int maxn = 5001, maxs = 5, maxp = 10;
int m, len, mod;
LL n;
struct Matrix
{
int r, c, num[maxs][maxs];
Matrix operator + (const Matrix &t) const
{
Matrix ret = {r, c};
for(int i = 0; i < r; ++i)
for(int j = 0; j < c; ++j)
if((ret.num[i][j] = num[i][j] + t.num[i][j]) >= mod)
ret.num[i][j] -= mod;
return ret;
}
Matrix operator * (const Matrix &t) const
{
Matrix ret = {r, t.c};
for(int i = 0; i < r; ++i)
for(int j = 0; j < c; ++j)
for(int k = 0; k < t.c; ++k)
ret.num[i][k] = (ret.num[i][k] + (LL)num[i][j] * t.num[j][k]) % mod;
return ret;
}
Matrix operator * (const int &t) const
{
Matrix ret = {r, c};
for(int i = 0; i < r; ++i)
for(int j = 0; j < c; ++j)
ret.num[i][j] = (LL)num[i][j] * t % mod;
return ret;
}
Matrix pow(LL k)
{
Matrix ret = {r, r}, tmp = *this;
for(int i = 0; i < r; ++i)
ret.num[i][i] = 1;
for( ; k; k >>= 1, tmp = tmp * tmp)
if(k & 1)
ret = ret * tmp;
return ret;
}
} I, A, B;
int mod_inv(int x)
{
return x <= 1 ? x : mod - mod / x * (LL)mod_inv(mod % x) % mod;
}
int mod_pow(int x, int k)
{
int ret = 1;
for( ; k; k >>= 1, x = (LL)x * x % mod)
if(k & 1)
ret = (LL)ret * x % mod;
return ret;
}
int origin()
{
int cnt = 0, tmp = mod - 1;
static int fact[maxp];
for(int i = 2; (LL)i * i <= tmp; ++i)
if(tmp % i == 0)
{
fact[cnt++] = i;
do tmp /= i; while(tmp % i == 0);
}
if(tmp > 1)
fact[cnt++] = tmp;
for(int ori = 2; ; ++ori)
{
bool flag = 1;
for(int i = 0; i < cnt && flag; ++i)
flag &= mod_pow(ori, (mod - 1) / fact[i]) != 1;
if(flag)
return ori;
}
}
int main()
{
for(int tot, sta, ter, w; scanf("%d%lld%d%d", &m, &n, &len, &mod) == 4; )
{
memset(&I, 0, sizeof(Matrix));
I.r = I.c = m;
for(int i = 0; i < m; ++i)
I.num[i][i] = 1;
memset(&A, 0, sizeof(Matrix));
A.r = A.c = m;
scanf("%d%d%d", &tot, &sta, &ter);
for(int u, v; tot--; )
{
scanf("%d%d", &u, &v);
++A.num[--u][--v];
}
memset(&B, 0, sizeof(Matrix));
B.r = B.c = m;
w = mod_pow(origin(), (mod - 1) / len);
for(int i = 0, ww = 1; i < len; ++i, ww = (LL)ww * w % mod)
B = B + (A * ww + I).pow(n);
printf("%d\n", (int)((LL)B.num[--sta][--ter] * mod_inv(len) % mod));
}
return 0;
}