BZOJ 2326: [HNOI2011]数学作业
自己手推一下然后矩阵就出来了。
可以发现对于矩阵:
\begin{pmatrix}
P_n\\
n\\
1
\end{pmatrix}
左乘矩阵(其中\( 10^k \)代表数位权):
\begin{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; }
然后。。。没了。。。
1 //{HEADS 2 #define FILE_IN_OUT 3 #define debug 4 #include <cstdio> 5 #include <cstring> 6 #include <cstdlib> 7 #include <cmath> 8 #include <ctime> 9 #include <algorithm> 10 #include <iostream> 11 #include <fstream> 12 #include <vector> 13 #include <stack> 14 #include <queue> 15 #include <deque> 16 #include <map> 17 #include <set> 18 #include <bitset> 19 #include <complex> 20 #include <string> 21 #define REP(i, j) for (int i = 1; i <= j; ++i) 22 #define REPI(i, j, k) for (int i = j; i <= k; ++i) 23 #define REPD(i, j) for (int i = j; 0 < i; --i) 24 #define STLR(i, con) for (int i = 0, sz = con.size(); i < sz; ++i) 25 #define STLRD(i, con) for (int i = con.size() - 1; 0 <= i; --i) 26 #define CLR(s) memset(s, 0, sizeof s) 27 #define SET(s, v) memset(s, v, sizeof s) 28 #define mp make_pair 29 #define pb push_back 30 #define PL(k, n) for (int i = 1; i <= n; ++i) { cout << k[i] << ' '; } cout << endl 31 #define PS(k) STLR(i, k) { cout << k[i] << ' '; } cout << endl 32 using namespace std; 33 #ifdef debug 34 #ifndef ONLINE_JUDGE 35 const int OUT_PUT_DEBUG_INFO = 1; 36 #endif 37 #endif 38 #ifdef ONLINE_JUDGE 39 const int OUT_PUT_DEBUG_INFO = 0; 40 #endif 41 #define DG if(OUT_PUT_DEBUG_INFO) 42 void FILE_INIT(string FILE_NAME) { 43 #ifdef FILE_IN_OUT 44 #ifndef ONLINE_JUDGE 45 freopen((FILE_NAME + ".in").c_str(), "r", stdin); 46 freopen((FILE_NAME + ".out").c_str(), "w", stdout); 47 48 #endif 49 #endif 50 } 51 typedef long long LL; 52 typedef double DB; 53 typedef pair<int, int> i_pair; 54 const int INF = 0x3f3f3f3f; 55 //} 56 57 LL n, p[20]; 58 /* { Matrix Multiplication index int型*/ 59 60 const int max_size = 4; 61 LL mod; 62 struct Matrix { 63 LL d[max_size][max_size]; 64 int r, c; 65 inline Matrix() { 66 memset(d, 0, sizeof d); 67 } 68 inline void set_size(int set_r, int set_c) { 69 r = set_r; 70 c = set_c; 71 } 72 inline void init(int size) { 73 r = c = size; 74 for(int i = 1; i <= r; ++i) { 75 d[i][i] = 1ull; 76 } 77 } 78 inline void print() { 79 for(int i = 1; i <= r; ++i) { 80 for(int j = 1; j < c; ++j) { 81 printf("%lld ", d[i][j]); 82 } 83 printf("%lld\n", d[i][c]); 84 } 85 } 86 }; 87 LL mul(LL a, LL b) { 88 LL ret = 0; 89 for(; b; b >>= 1, a = (a << 1) % mod) { 90 if(b & 1) { 91 ret = (ret + a) % mod; 92 } 93 } 94 return ret; 95 } 96 inline Matrix operator * (const Matrix &lhs, const Matrix &rhs) { 97 Matrix ret; 98 ret.set_size(lhs.r, rhs.c); 99 for(int i = 1; i <= lhs.r; ++i) { 100 for(int j = 1; j <= rhs.c; ++j) { 101 for(int k = 1; k <= lhs.c; ++k) { 102 ret.d[i][j] = (ret.d[i][j] + mul(lhs.d[i][k], rhs.d[k][j])) % mod; 103 } 104 } 105 } 106 return ret; 107 } 108 109 inline Matrix fast_pow(Matrix base, LL index) { 110 Matrix ret; 111 ret.init(base.r); 112 for(; index; index >>= 1, base = base * base) { 113 if(index & 1) { 114 ret = ret * base; 115 } 116 } 117 //ret.print(); 118 return ret; 119 } 120 121 122 /*} */ 123 124 int main() { 125 FILE_INIT("BZOJ2326"); 126 127 cin >> n >> mod; 128 p[0] = 1; 129 for(int i = 1; i <= 19; ++i) { 130 p[i] = (10 * p[i - 1]); 131 } 132 int len = 0; 133 for(LL tmp = n; tmp; tmp /= 10, ++len); 134 Matrix A, B; 135 A.set_size(3, 1); 136 A.d[1][1] = 0; 137 A.d[2][1] = 0; 138 A.d[3][1] = 1; 139 B.set_size(3, 3); 140 for(int i = 1; i < len; ++i) { 141 B.d[1][1] = p[i] % mod; B.d[1][2] = 1; B.d[1][3] = 1; 142 B.d[2][1] = 0; B.d[2][2] = 1; B.d[2][3] = 1; 143 B.d[3][1] = 0; B.d[3][2] = 0; B.d[3][3] = 1; 144 A = fast_pow(B, p[i - 1] * 9) * A; 145 // A.print(); 146 } 147 B.d[1][1] = p[len] % mod; B.d[1][2] = 1; B.d[1][3] = 1; 148 B.d[2][1] = 0; B.d[2][2] = 1; B.d[2][3] = 1; 149 B.d[3][1] = 0; B.d[3][2] = 0; B.d[3][3] = 1; 150 A = fast_pow(B, n - p[len - 1ull] + 1) * A; 151 cout << A.d[1][1] << endl; 152 153 return 0; 154 }
人的一切痛苦,本质上都是对自己的无能的愤怒。