BZOJ 2326: [HNOI2011]数学作业

自己手推一下然后矩阵就出来了。

可以发现对于矩阵：

egin{pmatrix}
P_n\
n\
1
end{pmatrix}

左乘矩阵（其中( 10^k )代表数位权）：

egin{pmatrix}
10^k & 1 &1 \
0& 1 & 1\
0& 0& 1
end{pmatrix}

即

 根据矩阵乘法的结合率，对于k分开考虑，进行快速幂计算即可。

有一个问题：能在计算时取模吗？不取模正确性是显然的，取模的话就比较dt了，不好说明答案没有改变。。如果有人知道如何证明取模是正确的，请回复sbit，万谢！

还有一个问题：两个大整数相乘中间结果可能爆掉，两个方法：一是把( 10^k )传入时先对m取模，二是模乘法，代码如下：

LL mul(LL a, LL b) {
    LL ret = 0;
    for(; b; b >>= 1, a = (a << 1) % mod) {
        if(b & 1) {
            ret = (ret + a) % mod;
        }
    }
    return ret;
}

然后。。。没了。。。

//{HEADS
#define FILE_IN_OUT
#define debug
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <algorithm>
#include <iostream>
#include <fstream>
#include <vector>
#include <stack>
#include <queue>
#include <deque>
#include <map>
#include <set>
#include <bitset>
#include <complex>
#include <string>
#define REP(i, j) for (int i = 1; i <= j; ++i)
#define REPI(i, j, k) for (int i = j; i <= k; ++i)
#define REPD(i, j) for (int i = j; 0 < i; --i)
#define STLR(i, con) for (int i = 0, sz = con.size(); i < sz; ++i)
#define STLRD(i, con) for (int i = con.size() - 1; 0 <= i; --i)
#define CLR(s) memset(s, 0, sizeof s)
#define SET(s, v) memset(s, v, sizeof s)
#define mp make_pair
#define pb push_back
#define PL(k, n) for (int i = 1; i <= n; ++i) { cout << k[i] << ' '; } cout << endl
#define PS(k) STLR(i, k) { cout << k[i] << ' '; } cout << endl
using namespace std;
#ifdef debug
#ifndef ONLINE_JUDGE
    const int OUT_PUT_DEBUG_INFO = 1;
#endif
#endif
#ifdef ONLINE_JUDGE
    const int OUT_PUT_DEBUG_INFO = 0;
#endif
#define DG if(OUT_PUT_DEBUG_INFO)
void FILE_INIT(string FILE_NAME) {
#ifdef FILE_IN_OUT
#ifndef ONLINE_JUDGE 
    freopen((FILE_NAME + ".in").c_str(), "r", stdin);
    freopen((FILE_NAME + ".out").c_str(), "w", stdout);
 
#endif
#endif
}
typedef long long LL;
typedef double DB;
typedef pair<int, int> i_pair;
const int INF = 0x3f3f3f3f;
//}
 
LL n, p[20];
/* { Matrix Multiplication index int型*/
 
const int max_size = 4;
LL mod;
struct Matrix {
    LL d[max_size][max_size];
    int r, c;
    inline Matrix() {
        memset(d, 0, sizeof d);
    }
    inline void set_size(int set_r, int set_c) {
        r = set_r;
        c = set_c;
    }
    inline void init(int size) {
        r = c = size;
        for(int i = 1; i <= r; ++i) {
            d[i][i] = 1ull;
        }
    }
    inline void print() {
        for(int i = 1; i <= r; ++i) {
            for(int j = 1; j < c; ++j) {
                printf("%lld ", d[i][j]);
            }
            printf("%lld
", d[i][c]);
        }
    }
};
LL mul(LL a, LL b) {
    LL ret = 0;
    for(; b; b >>= 1, a = (a << 1) % mod) {
        if(b & 1) {
            ret = (ret + a) % mod;
        }
    }
    return ret;
}
inline Matrix operator * (const Matrix &lhs, const Matrix &rhs) {
    Matrix ret;
    ret.set_size(lhs.r, rhs.c);
    for(int i = 1; i <= lhs.r; ++i) {
        for(int j = 1; j <= rhs.c; ++j) {
            for(int k = 1; k <= lhs.c; ++k) {
                ret.d[i][j] = (ret.d[i][j] + mul(lhs.d[i][k], rhs.d[k][j])) % mod;
            }
        }
    }
    return ret;
}
 
inline Matrix fast_pow(Matrix base, LL index) {
    Matrix ret;
    ret.init(base.r);
    for(; index; index >>= 1, base = base * base) {
        if(index & 1) {
            ret = ret * base;
        }
    }
    //ret.print();
    return ret;
}
     
 
/*} */
 
int main() {
    FILE_INIT("BZOJ2326");
 
    cin >> n >> mod;
    p[0] = 1;
    for(int i = 1; i <= 19; ++i) {
        p[i] = (10 * p[i - 1]);
    }
    int len = 0;
    for(LL tmp = n; tmp; tmp /= 10, ++len);
    Matrix A, B;
    A.set_size(3, 1);
    A.d[1][1] = 0;
    A.d[2][1] = 0;
    A.d[3][1] = 1;
    B.set_size(3, 3);
    for(int i = 1; i < len; ++i) {
        B.d[1][1] = p[i] % mod; B.d[1][2] = 1; B.d[1][3] = 1;
        B.d[2][1] =    0; B.d[2][2] = 1; B.d[2][3] = 1;
        B.d[3][1] =    0; B.d[3][2] = 0; B.d[3][3] = 1;
        A = fast_pow(B, p[i - 1] * 9) * A;
    //  A.print();
    }
    B.d[1][1] = p[len] % mod; B.d[1][2] = 1; B.d[1][3] = 1;
    B.d[2][1] =      0; B.d[2][2] = 1; B.d[2][3] = 1;
    B.d[3][1] =      0; B.d[3][2] = 0; B.d[3][3] = 1;
    A = fast_pow(B, n - p[len - 1ull] + 1) * A;
    cout << A.d[1][1] << endl;
 
    return 0;
}