• LightOJ 1132 Summing up Powers(矩阵快速幂)


    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

  • 相关阅读:
    Discuz 页面不能加载插件的原因和解决方法
    discuz 插件核心函数hookscript分析.
    比较容易犯的一些智障错误(不定时修改)
    浅谈树状数组入门
    图论的小总结
    usaco 2009 12 过路费
    0122(本来是想ak的但是因为智障只拿了200。)
    图论
    欧拉路
    bfs
  • 原文地址:https://www.cnblogs.com/LyonLys/p/LOJ_1132_Lyon.html
Copyright © 2020-2023  润新知