题解 P2020 【[NOI2011]兔农】

TimWYZ · 2020-08-14 17:01:15 · 题解

I 题意

给定整数 n, k, P

对于数组 F 有：

求 F_n \mod P

II 规律

假设 k = 7 . 我们写出 F 数组的前 30 项在模 k 意义下的值，每当出现 -1，就换行。

1 1 2 3 5 0
5 5 3 0
3 3 6 2 0
2 2 4 6 3 2 5 0 5 5 3 0
3 3 6 2 0
2 2 4 6 3 2 5 0 5 5 3 0
3 3 6 2 0
2 ...

可以发现，

3 3 6 2 0
2 2 4 6 3 2 5 0 5 5 3 0

就是一个循环节。在这个循环节当中，有几个性质

1. 对于循环节中的每一段，假设开头的数为 x ，那么数列为 x, x, 2x, 3x, 5x, 8x ... ，也是一个斐波那契数列。当遇到数列中某一项 \mod k = 1 时，数列结束
2. 从 1 可以得出，只要给定了数列开始的数，我们就能得知整个数列。
3. 一个段开始的数是上一个段的倒数第二个数。

通过这几个规律，我们就能搭建起这道题完整的解题框架了。

III 做法

既然有了循环节，我们就能通过循环节以及矩阵快速幂，快速算出 F_n\mod P。那么一个问题就是，如何寻求循环节？

在上面的规律中，我们发现，如果一个段开头是 x，长度为 len，普通的斐波那契数列的第 i 项为 fib_i ，那么这个段末尾的数（也就是要-1）的数就是 x * fib_{len}。我们要让它 x * fib_{len}\mod K = 1。看着是不是很眼熟？对，它就是乘法逆元。

因此，如果我们知道了 x，就能快速算出 fib_{len}，反之如果知道了 fib_{len}，也能快速算出 x。在这里我选择通过枚举 fib_{len} 计算 x。枚举斐波那契数列的每一项 fib_{i}，通过 exGCD 算出它的逆元 inv，对 k 取模得到 x，我们就能得出以 x 开头的段的长度为 i。

然后从 1 开始，每次计算以这个数为第一个数的段的末尾数，判断其之前是否出现过。如果出现，那么找到循环节；如果没出现，再将其作为段的开始迭代，直到找到循环节为止。

​ 特殊情况

​ 如果 a, b 不互质，那么 a 就没有 \mod b 意义下的逆元。在这道题目中，我们也没法保证每个 x 都有 \mod k 意义下的逆元。如果没有，那么情况将会如何呢？

​ 我们设 k = 10 （非质数），并写出它的前 30 项。

1 1 2 3 5 8 3 0
3 3 6 9 5 4 9 3 2 5 7 2 9 0
9 9 8 7 5 2 7 9 6 5 0
5 5 0 5 5 0 5 5 0 5 5 0 5 5 0 5 5       

​ 可以看出，在最末一行陷入了死循环，没有循环节。对于这种情况，我们要特殊判断，特殊处理。

得出每段长度 len 后，我们就可以着手推每个循环节的转移矩阵了。

对于一个行向量 mat_{i} 代表数列的第 i 项：

定义转移矩阵 tr1:

当需要 -1 时，定义转移矩阵 tr2:

那么，如果一个段的初始状态为 mat，长度为 len，那么它的结束状态就是 mat * tr1 ^ {len - 1} * tr2。 通过这一规律，我们就能算出每个循环节的总转移矩阵，然后计算出 F_n \mod P

注意：矩阵乘法全部是在 \mod P 意义下进行的。

另外一些细节详见代码。

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <iostream>
using namespace std;

typedef long long ll;
const int KR = 1e6 + 10, SZ = 3;

ll n, k, P;
ll kcnt = 0;
ll f[6 * KR], len[KR], seq[KR], vis[KR]; 
// f[i] 代表斐波那契第 i 位，len[i] 代表以 i 为开头的段的长度，seq 中存储所有段的开头，vis[i] 记录 i 作为开头的位置。
bool flag;

ll GCD(ll a, ll b) { // 由于还要判断是否存在逆元（互质），因此还需要一个朴素的GCD。
    if (!b) return a;
    return GCD(b, a % b);
}

void exGCD(ll a, ll b, ll &x, ll &y) {
    if (!b) {
        x = 1, y = 0;
        return;
    }
    exGCD(b, a % b, x, y);
    ll t = x;
    x = y;
    y = t - a / b * y;
}

ll getInv(ll a, ll P) {
    if (GCD(a, P) != 1) return -1; // 不互质，无逆元。
    ll x, y;
    exGCD(a, P, x, y);
    return (x % P + P) % P;
}

struct Matrix {
    ll o[SZ + 1][SZ + 1];

    Matrix() {
        memset(o, 0, sizeof(o));
    }

    Matrix operator * (const Matrix &x) const {
        Matrix ret;
        for (int i = 1; i <= SZ; i++)
            for (int j = 1; j <= SZ; j++)
                for (int k = 1; k <= SZ; k++)
                    ret.o[i][j] = (ret.o[i][j] + o[i][k] * x.o[k][j] + P) % P;
        return ret;
    }
} mat, tr1, tr2, tr;

Matrix quickPower(Matrix a, ll b) {
    Matrix ret;
    for (int i = 1; i <= SZ; i++) ret.o[i][i] = 1;

    while (b) {
        if (b & 1) ret = ret * a;
        a = a * a;
        b >>= 1;
    }
    return ret;
}

int main()
{
    freopen("luogu.in", "r", stdin);
    freopen("luogu.out", "w", stdout);

    scanf("%lld%lld%lld", &n, &k, &P);

    if (n == 1 || n == 2) {
        printf("1\n");
        return 0;
    }

    memset(len, 999999, sizeof(len));
    f[1] = f[2] = 1;

    for (ll i = 3; ; i++) {
        f[i] = (f[i - 1] + f[i - 2]) % k;
        if (f[i] % k == 1 && len[1] > 1e18) len[1] = i; // 需要特殊计算 1 开头的长度
        if (f[i] == 1 && f[i - 1] == 1) break;

        ll inv = getInv(f[i], k);
        if (inv != -1) len[inv % k] = min(len[inv % k], i);
    }

    ll now = 1, tot = 0;
    while (1) {
        seq[++kcnt] = now;
        vis[now] = kcnt;
        if (len[now] > 1e18) { // 如果没有逆元
            for (int i = 1; i < kcnt; i++) tot += len[seq[i]]; // 计算之前段的总长
            flag = 1;
            break;
        }

        now = (now * f[len[now] - 1]) % k;
        if (vis[now]) { // 这个数之前出现过，发现循环节
            for (int i = 1; i < vis[now]; i++) tot += len[seq[i]]; // 计算循环节之前段的总长
            break;
        }
    }

    mat.o[1][1] = mat.o[1][3] = 1;
    tr1.o[1][1] = tr1.o[1][2] = tr1.o[2][1] = tr1.o[3][3] = 1;
    tr2.o[1][1] = tr2.o[1][2] = tr2.o[2][1] = tr2.o[3][3] = 1;
    tr2.o[3][1] = -1;

    if (n <= tot) { // 要特别判断如果 n 在无逆元序列 / 循环节之前的情况（两个样例都是这种情况），暴力一段一段算
        len[1]--; // 初始已经算出了 F[1]，因此第一段的长度要 -1 ，n 随之也 -1
        n--;
        for (int i = 1; i < vis[now]; i++) {
            if (n >= len[seq[i]]) { // 假设 n 还够完整的一段，直接矩阵快速幂
                mat = mat * quickPower(tr1, len[seq[i]] - 1) * tr2; 
                n -= len[seq[i]];
            } else { // 否则暴力乘
                mat = mat * quickPower(tr1, n);
                printf("%lld\n", mat.o[1][1]);
                return 0;
            }
        }
    } else {
        len[1]--;
        n -= tot;

        for (int i = 1; i < vis[now]; i++) mat = mat * quickPower(tr1, len[seq[i]] - 1) * tr2; // 此处要计算出进入循环节前的状态

        if (flag) { // 无逆元情况，直接用 tr1 快速幂
            mat = mat * quickPower(tr1, n);
            printf("%lld\n", mat.o[1][1]);
        } else {
            ll loopLen = 0;
            for (ll i = 1; i <= SZ; i++) tr.o[i][i] = 1;
            for (ll i = vis[now]; i <= kcnt; i++) { // 否则计算循环节的总转移矩阵 tr
                tr = tr * quickPower(tr1, len[seq[i]] - 1) * tr2;
                loopLen += len[seq[i]];
            }
            ll tmp = n / loopLen; // 算一下 n 中有几个完整的循环节
            mat = mat * quickPower(tr, tmp);
            n = n - loopLen * tmp; // 剩下的部分暴力算

            for (ll i = vis[now]; i <= kcnt; i++) {
                if (n >= len[seq[i]]) { // 假设 n 还够完整的一段，直接矩阵快速幂
                    mat = mat * quickPower(tr1, len[seq[i]] - 1) * tr2;
                    n -= len[seq[i]];
                } else { // 否则暴力乘
                    mat = mat * quickPower(tr1, n);
                    printf("%lld\n", mat.o[1][1]);
                    return 0;
                }
            }
        }

    }

    return 0; 
}