Jan's LightOJ :: Problem 1132 - Summing up Powers
我的队友xdp给我做的一道矩阵快速幂的题。昨晚花了两三个小时,想出了怎么用矩阵来实现转移,今天早上校赛前打了预处理,然后校赛回来以后就把它A了!
题意就是求给出式子的结果mod 2^32的最终结果。
开始的时候,还真的想不到怎么用矩阵实现快速的状态转移,只好百度一下这个式子有没有什么公式。不过,一轮搜索过后还是没有可以套矩阵的方法,或者可以快速求解的公式,所以只好把这个式子写成关于n的递推式,然后自己观察这个式子每次的变化。在观察的时候,我尝试着每两个(n-1)^k和n^k求差,如果得到的不是常数列,就继续将得到的序列相邻元素两两相减,直到找到这个常数列为止。找到以后,发现最终的常数列都是k!,恍然大悟,这个跟求导的做法很相似。当时我还认真的研究了一下,什么时候可以利用矩阵,其实也就是要满足转移条件的情况下,要能有一个统一的转移方法。
借助不断两两相减,得到常数列,这个过程的每一个序列的第一个数都取出来就是我们要的初始序列了。然后,逆着求差的思路,借助矩阵,我们可以得到一个递推矩阵。接着的矩阵快速幂就简单了。另外,因为求模的是2^32,所以只需要直接的用unsigned int,然后任由它溢出就好了。代码如下:
View Code
1 #include <cstdio> 2 #include <cstring> 3 #include <algorithm> 4 #include <iostream> 5 6 using namespace std; 7 8 typedef unsigned int UI; 9 typedef long long LL; 10 const int N = 55; 11 UI coi[N][N][N]; 12 13 #define _clr(x) memset(x, 0, sizeof(x)) 14 #define REP_1(i, n) for (int i = 1; i <= (n); i++) 15 #define REP(i, n) for (int i = 0; i < (n); i++) 16 17 void PRE() { 18 _clr(coi); 19 REP_1(i, 53) { 20 REP(j, i + 1) { 21 coi[i][0][j] = (UI) j; 22 REP(k, i - 1) { 23 coi[i][0][j] *= (UI) j; 24 } 25 } 26 REP_1(j, i) { 27 REP(k, i) { 28 coi[i][j][k] = coi[i][j - 1][k + 1] - coi[i][j - 1][k]; 29 } 30 } 31 // REP(j, i + 1) { 32 // REP(k, i + 1) cout << coi[i][j][k] << ' '; 33 // cout << endl; 34 // } 35 } 36 } 37 38 const int maxSize = 55; 39 int curSize = maxSize; 40 41 struct Mat { 42 UI val[maxSize][maxSize]; 43 Mat(UI one = 0) { 44 REP(i, curSize) { 45 REP(j, curSize) { 46 val[i][j] = 0; 47 } 48 val[i][i] = one; 49 } 50 } 51 void print() { 52 cout << "mat" << endl; 53 REP(i, curSize) { 54 REP(j, curSize) { 55 cout << val[i][j] << ' '; 56 } 57 cout << endl; 58 } 59 cout << "!!!" << endl; 60 } 61 } Base, op; 62 63 Mat operator * (Mat &_a, Mat &_b) { 64 Mat ret = Mat(); 65 REP(i, curSize) { 66 REP(k, curSize) { 67 if (_a.val[i][k]) { 68 REP(j, curSize) { 69 ret.val[i][j] += _a.val[i][k] * _b.val[k][j]; 70 } 71 } 72 } 73 } 74 return ret; 75 } 76 77 Mat operator ^ (Mat &__a, LL _p) { 78 Mat _a = __a; 79 Mat ret = Mat(1); 80 while (_p > 0) { 81 if (_p & 1) ret = ret * _a; 82 _a = _a * _a; 83 _p >>= 1; 84 } 85 return ret; 86 } 87 88 void makeMat(int n) { 89 Base = Mat(); 90 op = Mat(); 91 curSize = n + 2; 92 REP_1(i, n) { 93 Base.val[0][i] = coi[n][i - 1][0] + coi[n][i][0]; 94 } 95 Base.val[0][n + 1] = coi[n][n][0]; 96 // Base.print(); 97 REP(i, curSize) { 98 REP(j, 2) { 99 op.val[i + j][i] = 1; 100 } 101 } 102 // op.print(); 103 } 104 105 UI work(LL n, int k) { 106 makeMat(k); 107 // Base.print(); 108 op = op ^ n; 109 Base = Base * op; 110 return Base.val[0][0]; 111 } 112 113 const LL mask = 0xffffffffll; 114 115 UI bruceForce(LL n, int k) { 116 UI ret = 0; 117 REP_1(i, n) { 118 UI tmp = 1; 119 REP(j, k) { 120 tmp *= (UI) i; 121 } 122 ret += tmp; 123 } 124 return ret; 125 } 126 127 int main() { 128 // freopen("in", "r", stdin); 129 PRE(); 130 int T, k; 131 LL n; 132 cin >> T; 133 REP_1(cas, T) { 134 cin >> n >> k; 135 printf("Case %d: ", cas); 136 if (k) cout << work(n, k) << endl; 137 else cout << (n & mask) << endl; 138 // cout << bruceForce(n, k) << endl; 139 } 140 return 0; 141 }
——written by Lyon