OpenJudge 1058 Guideposts | BZOJ 3328 PYXFIB

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

题意:
两题的题面不太相同,抽象后题意如下。
给出一个M×M的矩阵G,求

i=0N[K|i](Ni)Gis,t(modP)

其中P是质数,并且P1(modK)
Guideposts:1N1018,1M5,1K5000,1P109
PYXFIB:1N1018,M=2,1K2104,1P109

题解:
不考虑[K|i]的话,有Ni=0(Ni)Gis,t=Ni=0(Ni)(INiGi)s,t=(I+G)Ns,t
[K|i]可以利用单位复根进行构造,显然有

0j<KωijKK=[K|i]

则原式可化为
i=0Nj=0K1ωijKK(Ni)Gis,t=1Kj=0K1i=0N(Ni)(ωjKG)is,t=1Kj=0K1i=0N(Ni)(INi(ωjKG)i)s,t=1Kj=0K1(I+ωjKG)Ns,t

因此可以得到一个O(M3KlogN)的做法,但是直接使用单位复根,会导致精度过低。
由于P是质数,并且P1(modK),即K|φ(P),故存在x使得ord(x)=K,那么ωKx(modP),因此可以找出一个原根g,可以得到ord(gφ(P)K)=K,于是用gφ(P)K代替ω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;
}
©️2020 CSDN 皮肤主题: 技术黑板 设计师: CSDN官方博客 返回首页
实付0元
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值