OpenJudge 1058 Guideposts | BZOJ 3328 PYXFIB

题目：
http://poj.openjudge.cn/practice/1058/
http://www.lydsy.com/JudgeOnline/problem.php?id=3328

题意：
两题的题面不太相同，抽象后题意如下。
给出一个$M\times M$的矩阵$G$，求

$$\sum_{i=0}^{N}{[K|i]\binom{N}{i}G^i_{s,t}}\pmod{P}$$
其中$P$是质数，并且$P\equiv 1\pmod{K}$。
Guideposts：$1\leq N\leq 10^{18},1\leq M\leq 5,1\leq K\leq 5000,1\leq P\leq 10^9$
PYXFIB：$1\leq N\leq 10^{18},M=2,1\leq K\leq 2\cdot 10^4,1\leq P\leq 10^9$

题解：
不考虑$[K|i]$的话，有$\sum_{i=0}^{N}{\binom{N}{i}G^i_{s,t}}=\sum_{i=0}^{N}{\binom{N}{i}(I^{N-i}G^i)_{s,t}}=(I+G)^N_{s,t}$。
而$[K|i]$可以利用单位复根进行构造，显然有

$$\frac{\sum_{0\leq j<K}{\omega_K^{ij}}}{K}=[K|i]$$
则原式可化为 $$\sum_{i=0}^{N}\sum_{j=0}^{K-1}{\frac{\omega_K^{ij}}{K}\binom{N}{i}G^i_{s,t}}\\=\frac{1}{K}\sum_{j=0}^{K-1}\sum_{i=0}^{N}{\binom{N}{i}(\omega_K^{j}G)^i_{s,t}}\\=\frac{1}{K}\sum_{j=0}^{K-1}\sum_{i=0}^{N}{\binom{N}{i}(I^{N-i}(\omega_K^{j}G)^i)_{s,t}}\\=\frac{1}{K}\sum_{j=0}^{K-1}{(I+\omega_K^{j}G)^N_{s,t}}$$
因此可以得到一个 $O(M^3K\log N)$的做法，但是直接使用单位复根，会导致精度过低。
由于 $P$是质数，并且$P\equiv 1\pmod{K}$，即 $K|\varphi(P)$，故存在 $x$使得$ord(x)=K$，那么 $\omega_K\equiv x\pmod{P}$，因此可以找出一个原根 $g$，可以得到$ord(g^\frac{\varphi(P)}{K})=K$，于是用 $g^\frac{\varphi(P)}{K}$代替 $\omega_K$即可。

代码：

#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;
}