Matrix Power Series

Posted by twcai in Geek Stuff on May 25, 2009

经小HH大牛指点，重新过了一遍POJ 3233，效率有很大的提高
题目给出一个矩阵A, 要求求出A^1 + A^2 + … + A^k
这里的做法是构造一个两倍于A的矩阵Bmatrix{2}{2}{e A 0 A}
然后求B^K次即可，所要的矩阵位于B的右上角

  1. #include <iostream>
  2. using namespace std;
  3. //A^1 + A^2 + ... + A^k
  4. int modnum;
  5.  
  6. class matrix{
  7. public :
  8.     int rw, cn, ma[65][65];
  9.  
  10.     matrix() {
  11.         rw = 0, cn = 0;
  12.     }
  13.  
  14.     matrix(int x, int y) {
  15.         int i, j;
  16.         rw = x;
  17.         cn = y;
  18.         for(i = 0; i < rw; i++) {
  19.             for(j = 0; j < cn; j++) { ma[i][j] = 0; }
  20.         }
  21.     }
  22.  
  23.     matrix(const matrix &x) {
  24.         int i, j;
  25.         rw = x.rw, cn = x.cn;
  26.         for(i = 0; i < rw; i++){
  27.             for(j = 0; j < cn; j++) { ma[i][j] = x.ma[i][j]; }
  28.         }
  29.     }
  30.  
  31.     matrix operator*(matrix &x)
  32.     {
  33.         int i, j, k;
  34.         matrix ret(cn, x.rw);
  35.         for(i = 0; i < rw ; i++) {
  36.             for(j = 0; j < x.cn ; j++) {
  37.                 for(k = 0; k < cn; k++) {
  38.                     ret.ma[i][j] += ma[i][k] * x.ma[k][j];
  39.                     if(ret.ma[i][j] >= modnum)
  40.                         ret.ma[i][j] %= modnum;
  41.                 }
  42.             }
  43.         }
  44.         return ret;
  45.     }
  46.  
  47.     matrix power(int x) {
  48.         matrix ret(rw, cn);
  49.         if(x == 1) return *this;
  50.         if(x == 2) {
  51.             int i, j ,k;
  52.             for(i = 0; i < rw; i++) {
  53.                 for(j = 0; j < cn; j++) {
  54.                     for(k = 0; k < cn; k++) {
  55.                         ret.ma[i][j] += ma[i][k] * ma[k][j];
  56.                         if(ret.ma[i][j] >= modnum)
  57.                             ret.ma[i][j] %= modnum;
  58.                     }
  59.                 }
  60.             }
  61.         }
  62.         else if(!(x & 1)) {
  63.             ret = this->power(x >> 1);
  64.             ret = ret.power(2);
  65.         }
  66.         else {
  67.             ret = this->power(x >> 1);
  68.             ret = ret.power(2);
  69.             ret = ret * (*this);
  70.         }
  71.         return ret;
  72.     }
  73.  
  74.     inline matrix operator+= (const matrix &x) {
  75.         int i, j;
  76.         for(i = 0; i < rw; i++) {
  77.             for(j = 0; j < cn; j++)
  78.             {
  79.                 ma[i][j] += x.ma[i][j];
  80.                 if(ma[i][j] >= modnum)
  81.                     ma[i][j] %= modnum;
  82.             }
  83.         }
  84.         return *this;
  85.     }
  86.  
  87.     inline matrix operator = (const matrix &x) {
  88.         int i, j;
  89.         rw = x.rw, cn = x.cn;
  90.         for(i = 0; i < rw; i++) {
  91.             for(j = 0; j < cn; j++) { ma[i][j] = x.ma[i][j]; }
  92.         }
  93.         return *this;
  94.     }
  95. };
  96.  
  97. int main()
  98. {
  99.     int n, m, k;
  100.     int i, j;
  101.     while(scanf("%d %d %d", &n, &k, &m) != EOF) {
  102.         matrix org(2 * n, 2 * n), ret(2 * n, 2 * n);
  103.         for(i = 0; i < n; i++) {
  104.             org.ma[i][i] = 1;
  105.             for(j = 0; j < n; j++) {
  106.                 scanf("%d", &org.ma[i][n + j]);
  107.                 org.ma[n + i][n + j] = org.ma[i][n + j];
  108.             }
  109.         }
  110.         modnum = m;
  111.         ret = org.power(k);
  112.         for(i = 0; i < n; i++) {
  113.             for(j = 0; j < n; j++) {
  114.                 if(j) printf(" ");
  115.                 printf("%d", ret.ma[i][n + j]);
  116.             }
  117.             printf("\n");
  118.         }
  119.     }
  120.     return 0;
  121. }

2 Comments

省赛总结

Posted by twcai in ACM/ICPC on May 24, 2009

rank name solved A B C D E F G H I J K 罚时
15 HDU-Firebirds 7 13 (1) 144 (3) 138 (1) 7 0 22 (1) 0 0 24 (1) 121 (2) 66 (2) 608

省赛对我来说是入队以来最重要的比赛，因为直到那天我都还没想好省赛后是否退役，说不定这就是我最后一场比赛了(明年省赛不算，如果退役了，相信那时候也只是抱着旅游的态度了)。省赛前夜11点就整理完东西，洗漱上床了，但一直醒到2点多。

第二天6点半起来了，感觉有点睡眠不足，不过还算蛮精神的。上午的东西都直接略过，下午在机房门口等比赛开始的时候，发现去年杭州赛区的主场作战真的是蛮大的优势。

比赛开始我猜C是ms题，然后看了题目，是跟通常情况有点不同的MST，我们队不熟prim，又不会证明kruscal肯定可以得到正解，于是我放下C看I，并发现是简单题上机敲。敲完sample出不来，于是打印。Teddy上机敲A，在无声无息中过掉了。我发现自己I想复杂了，上去改了下，还是有问题。甘露上机敲F，期间我发现I的一个bug，甘露过掉F后，我上去改了I，提交后Yes。

过完I，看了下statistics，B和K都很多人过，我推了下B，是一个B - (X + A / X)的函数，忘掉怎么求最小值，直接扔给斌哥了。然后看K，简单构造题，上机打完后，Teddy对我的解法不是很肯定，造了几组数据测了下，然后提交CE……崩溃，赛前测编译器，好像就是平常poj的g++编译器，没想到还是出现这种错误。然后加了头文件，返回Yes。wyb上机打B。

接着看了下J，很熟悉的题目，记得是一个DP。不过Teddy对于这题的解法记忆不是很清晰了，我想了个简单的解法，换下wyb实现了下，过不了造的数据。这时teddy记起来解法了，换下wyb敲J。我帮斌哥一起看B，另外想了个解法，不过到最后还是跟斌哥已经wa的做法一样，感觉有两种可能，一种是情况没考虑全，一种是精度问题，然后一直在讨论。teddy过了J后，对C还是拿不定，我觉得还是应该试下kruscal。甘露很快敲好了C，提交以后居然Yes了……昏倒。接着wyb上去改了一下B，也过掉了。

这时候过了两个钟头多，我们的rank好像是在11，但是罚时没有优势。我们都很清楚要稳拿金牌一定要再过一题，而且剩下时间还蛮充裕的。wyb最先开了D题，我和甘露讨论了G。但是我图论实在烂，在读过几乎全部题目后(除了D)决定用搜索开H。wyb用暴力打好D，提交后tle，然后甘露上机敲G，敲到一半他觉得不肯定，换我上机敲H  =,=||

赛后想起来这段时间真的是在梦游，我敲H的时候心里很虚，但是还是觉得不应该放弃有想法的题目，敲完大部分代码。敲到没想清楚的地方后，换上wyb继续敲D。我记得这时候已经剩下一个钟头不到的时间了，于是和teddy都渐渐放弃了自己的题目，和wyb一起做D。这题他用了上次宁波刚从DD牛那里学到的逆元来处理做除法取模，本来我们都以为这题会很有希望的，没想到就一直这么修修改改到比赛结束。

时间还是无情的指向了5:25。我们已经明白金牌无望，伤心下一直在抱怨题目区分度不好。但是冷静下来想想，这真的只是题目的原因吗？我们队还是在比赛策略上有很大的失误：

一是C题太犹豫。Kruscal只需要5分钟的机时，但是我们一直到2个多钟头才下决心去试，是对比赛形势的错误估计。这种比赛，有6、7个队过掉了就应该果断的尝试。

二是7题以后的各自为战。我们队的优势是互补，平常比赛难题也都是合作解出的。

三是最后合力攻D时，阵脚已经乱了。两个人不去好好读题，却不停修改代码以及问wyb各种问题，影响了他深入思考。到比赛结束teddy还没有理解题意，我还没有读懂wyb的解法，最后一次提交代码的版本都搞乱了。这种时候，应该有一个人出来稳定整个队的情绪。

四是练习赛没有把握好测环境的机会，导致CE的出现。这是经验和总结的欠缺。

银牌实在是令人遗憾，但也是真实水平的反映。每一次比赛都是绝好的积累经验和反省的机会，希望再通过一个暑假的集训，Firebirds能拿到让我们三个人都满意的成绩。

No Comments

代码量记录

Posted by twcai in Geek Stuff on May 23, 2009

至今为止我保存下来的代码一共有
Update at 2008-07-05 : 13547 行 256347 字节
Update at 2008-07-13 : 14749 行 278236 字节
Update at 2008-07-26 : 15858 行 301407 字节
Update at 2008-08-03 : 16605 行 312928 字节
Update at 2008-08-12 : 17728 行 332990 字节
Update at 2008-08-27 : 19169 行 360280 字节 (这段时间很松懈，代码量很少，要加把劲才行)
Update at 2008-09-05 : 20638 行 388631 字节
Update at 2008-10-03 : 22766 行 427148 字节 (不行啊我，开学了还是要抓紧)
Update at 2008-10-21 : 24384 行 456496 字节 (+U+U)
Update at 2008-10-28 : 26128 行 483370 字节 (这段时间废话太多)
Update at 2008-11-04 : 29260 行 536219 字节 (加上了一些以前没写完的代码)
Update at 2008-11-13 : 30342 行 553013 字节 (Regional 加油！)
Update at 2008-11-26 : 32468 行 592289 字节 (要去成都了，淡定～)
Update at 2008-12-21 : 33901 行 612669 字节 (期末好忙啊，时间本来就少还都被我浪费了)
Update at 2009-02-05 : 35432 行 634096 字节 (一个慵懒的寒假)
Update at 2009-02-23 : 38100 行 684674 字节
Update at 2009-03-01 : 40051 行 725150 字节 (代码是在写着，但是很多都只完成了一半)
Update at 2009-03-07 : 41218 行 745806 字节 (离上次刚好一个星期)
Update at 2009-03-18 : 43855 行 798499 字节 (选拔赛和省赛近在眼前了，加油！)
Update at 2009-03-25 : 45081 行 816814 字节 (最近没做什么有质量的题目，这些代码基本来自水题)
Update at 2009-04-05 : 46283 行 838456 字节 (太懒了，一定要严格要求！)
Update at 2009-04-17 : 47278 行 853792 字节 (越来越锉了……)
Update at 2009-05-07 : 48699 行 890131 字节 (省赛！！一直要完成专题)
Update at 2009-05-22 : 50731 行 931529 字节 (终于5万行。参加比赛一年了，明天又是去ZJG的日子，加油啊)

这篇日志从去年暑假集训刚开始开始记录，一直设为私有。这一年来的生活算是落幕了，接下去是新的开始。

2 Comments

省赛前夕

Posted by twcai in ACM/ICPC on May 18, 2009

省赛就要到了，组队选拔赛也终于在上个周末全部结束。
最近的四场比赛我们队的发挥没有前段时间好了，不过还是有进步的地方，趁省赛前回头看看，总结一下。
可能是之前的排名给了我们压力，这几场周天涯队的发挥回复正常后，明显我们在开局变得有些急躁了。我感觉自己进入比赛状态比较慢，其他队AC会分散我的注意力，导致这几场简单题都出得慢。到了省赛上，牛队更多，自己的心态还需要调节吧。不过有一件值得乐观的事情就是我们队打破了比赛最后一小时不出题的坏习惯，我相信这会给我们队在正式比赛中带来信心。
开学到现在就一直在看数论，进展很慢，不过还是感到有点进展了，可能数论的知识结构的建立本来就需要时间吧，再加上我学习习惯这么烂……汗个。省赛前已经没时间看了，要回过头复习下以前的知识点。
这个学期就一直在犹豫省赛后还要不要继续。周末阿光来看老婆，于是和他交流了下，然后这两天跟Teddy和刘老师都交流了下，我越来越偏向继续了。不过决定我会在跟浙大的老师请教过后再下，这样也算给家里一个交代。
恩，我越来越期待省赛后的日子了～

No Comments

数论专题

Posted by twcai in ACM/ICPC on May 18, 2009

要开数论专题了，发一篇日志报告做题进度来督促自己。
知识点：
独立剩余

进位制： 进位制非常强大，在两个64位整数之间运算并取模的过程中，可以防止溢出(不需要动用大数)，而且复杂度比较低。
来看乘法运算和幂运算。设x, y , n是三个long long 的整数，我们需要分别计算 (x * y) % n 和 (x ^ y) % n .
首先是乘法，把 y 转化成二进制，y = sum{i = 0}{k}{a_i * 2 ^ i} , 那么 x * y = sum{i = 0}{k}{(x * a_i * 2 ^ i) % n} . 因为a_i的取值只有0 和 1, 所以这个计算的过程被大大简化了。
幂运算基本做法相同，x ^ y = prod{i = 0}{k}{(x * a_i * 2 ^ i) % n} , 在幂运算计算过程中 x * a_i * 2 ^ i 可以由前一个推得，并且需要使用到乘法运算。

Update : 本来看数论就是因为受了宁波理工A题的刺激，今天在黄队的读书笔记上终于找到这个知识点 “独立剩余”

POJ:
1006
5099752 clumsydragon 1006 Accepted 392K 32MS G++ 844B
中国剩余定理
1061
5082504 clumsydragon 1061 Accepted 400K 16MS G++ 870B
扩展欧几里德算法解线性方程 : A * (m - n) + V * L = y - x (m > n)
结果就是大于0的第一个解 A_i * (y - x) / {gcd(m - n, L)}
1405
1423
5090074 clumsydragon 1423 Accepted 8236K 282MS G++ 699B
使用公式 n! = sqrt{2 * n^2} * (n/e)^n * (1 + 1/(12*n) + 1/(288* n^2) + O(1/n^3))
据说是《计算机程序设计艺术》中的一个公式，好强大。
1606
1811
Miller_Rabin + 快速幂取模 + Pollard的rho启发式整数因子分解
1845
2447
2700
3146
3219
3358

HDOJ:
1930
1348734 2009-05-08 13:57:06 Accepted 1930 0MS 268K 1182 B G++ CD@HDU
比较单纯的中国剩余定理

No Comments