洛谷P3216 [HNOI2011]数学作业 题解 矩阵快速幂

题目链接：https://www.luogu.com.cn/problem/P3216

题目大意：定义函数 \(f(n)\) 为将 \(1 \sim n\) 连起来得到的数，求 \(f(n) \text{ mod } m\) 的结果。

解题思路：

状态转移方程为：

\[f_i = 10^{ \lfloor \log_{10} i \rfloor } \times f_{i-1} + i \]

转移矩阵为：

\[\begin{bmatrix} f_n \\ n \\ 1 \end{bmatrix} = \begin{bmatrix} 10^k & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 1 \end{bmatrix} \times \begin{bmatrix} f_{n-1} \\ n-1 \\ 1 \end{bmatrix} \]

其中的 \(k\) 表示 \(\lfloor \log_{10} i \rfloor + 1\)（即 \(n\) 所占的位数）。但是 \(k\) 是变化的，没有办法直接矩阵乘法。但是可以发现：

• 当 \(n \in [1, 9]\) 时，\(k = 1\)；
• 当 \(n \in [10, 99]\) 时，\(k = 2\)；
• 当 \(n \in [100, 999]\) 时，\(k=3\)；
• …… ……
• 当 \(n = 10^{18}\) 时，\(k=19\)。

所以 \(k\) 总共有 \(19\) 种不同的情况，我们只需要针对不同的 \(k\) 做矩阵快速幂就可以了。

示例代码：

#include <bits/stdc++.h>
using namespace std;
long long n, m;
struct Matrix {
    long long a[3][3];
    Matrix operator * (Matrix b) const {
        Matrix c;
        memset(c.a, 0, sizeof(c.a));
        for (int i = 0; i < 3; i ++)
            for (int j = 0; j < 3; j ++)
                for (int k = 0; k < 3; k ++)
                    c.a[i][j] = (c.a[i][j] + a[i][k] * b.a[k][j] % m) % m;
        return c;
    }
    Matrix operator ^ (long long n) const {
        Matrix b, c;
        memcpy(b.a, a, sizeof(a));
        memset(c.a, 0, sizeof(c.a));
        for (int i = 0; i < 3; i ++) c.a[i][i] = 1;
        while (n) {
            if (n % 2) c = c * b;
            b = b * b;
            n /= 2;
        }
        return c;
    }
} ans, tmp;
long long x[20];
int main() {
    cin >> n >> m;
    long long t1 = 1, t2 = 9;
    for (int i = 0; i < 18; i ++) {
        x[i] = min(t2, n) - t1 + 1;
        if (n < t2) break;
        t1 *= 10, t2 = t2 * 10 + 9;
    }
    if (n == 1000000000000000000LL) x[18] = 1;
    memset(ans.a, 0, sizeof(ans.a));
    for (int i = 0; i < 3; i ++) ans.a[i][i] = 1;
    for (int i = 18; i >= 0; i --) {
        if (!x[i]) continue;
        memset(tmp.a, 0, sizeof(tmp.a));
        tmp.a[0][0] = 1;
        for (int j = 0; j <= i; j ++) tmp.a[0][0] = tmp.a[0][0] * 10 % m;
        tmp.a[0][1] = tmp.a[0][2] = tmp.a[1][1] = tmp.a[1][2] = tmp.a[2][2] = 1;
        ans = ans * (tmp ^ x[i]);
    }
    cout << ans.a[0][2] << endl;
    return 0;
}