熟练掌握动态规划

---

状态压缩DP

---

旅行商问题

• 问题描述：给定一个n个顶点组成的带权有向图的距离矩阵d(i,j)，要求从顶点0出发，经过每个顶点恰好一次后再回到顶点0.问所经过的边的总权重的最小值是多少？
• 限制条件：
  2≤n≤15
  0≤d(i,j)≤1000
• 分析：著名的旅行商问题(TSP, Traveling Salesman Problem)。TSP问题是NP困难的，没有已知的多项式时间的高效算法。
  所有可能的路线共有(n-1)!种，无法试遍每一种情况。
  对于这个问题，我们可以使用DP。先写递推式，假设现在已经访问过的顶点的集合(起点0当作还未访问过的顶点)为S，当前所在的顶点为v，用dp[S][v]表示从v出发访问剩余的所有顶点，最终回到顶点0的路径的权重总和的最小值。由于从v出发可以移动到任意的一个节点u∉S，因此有如下递推式：
  dp[V][0]=0
  dp[S][v]=min{dp[ S∪{u} ] [u] + d(v,u) |u∉S}
  我们只要按照这个递推式进行计算就可以了，由于在这个递推式中有一个下标是集合而不是普通的整数，因此需要稍加处理。首先我们试着使用记忆化搜索求解，虽然有一个下标不是整数，但我们可以把它编码成一个整数，或者个它定义一个全序关系并用二叉搜索树存储，从而可以使用记忆化搜索来求解。特别地，对于集合我们可以把每一个元素的选取与否对应到一个二进制位里，从而把状态压缩成一个整数，大大方便了计算和维护。
• 代码： 

  #include <cstdio>
  #include <cstdio>
  #include <cctype>
  #include <cstring>
  #define number s-'0'
  
  using namespace std;
  
  const int MAX_N=20;
  const int INF=0x3f3f3f3f;
  int n,m;
  int d[MAX_N][MAX_N];
  int dp[1<<MAX_N][MAX_N];
  int pre[MAX_N];
  
  void read(int &x)
  {
      char s;
      x=0;
      bool flag=0;
      while(!isdigit(s=getchar()))
          (s=='-')&&(flag=true);
      for(x=number;isdigit(s=getchar());x=x*10+number);
      (flag)&&(x=-x);
  }
  
  void write(int x)
  {
      if(x<0)
      {
          putchar('-');
          x=-x;
      }
      if(x>9) write(x/10);
      putchar(x%10+'0');
  }
  
  int rec(int S, int v);
  
  int main()
  {
      read(n);read(m);
      memset(d, 0x3f, sizeof(d));
      for (int i=0; i<n; i++) i[i[d]]=0;
      for (int i=0; i<m; i++)
      {
          int x,y,z;
          read(x);read(y);read(z);
          d[x][y]=z;
      }
      memset(dp, -1, sizeof(dp));
      write(rec(0,0));
      putchar('
  ');
  }
  
  int min(int x, int y)
  {
      if (x<=y) return x;
      return y;
  }
  
  int rec(int S, int v)
  {
      if (dp[S][v]>=0) return dp[S][v];
      if (S==(1<<n)-1 && v==0) return dp[S][v]=0;
      int res=INF;
      for (int u=0; u<n; u++)
      {
          if (!(S>>u&1)) res=min(res, rec(S|1<<u, u)+d[v][u]);
      }
      return dp[S][v]=res;
  }

   这样，就可以在O(2ⁿn²)的时间内完成计算。对于不是整数的情况，很多时候很难确定一个合适的递推顺序，因此使用记忆化搜索可以避免这个问题。
  不过，在这个问题中，对于任意两个整数i和j，如果它们对应的集合满足S(i)⊆S(j)，就有i≤j，因此还可通过循环求出答案：

  #include <cstdio>
  #include <cstdio>
  #include <cctype>
  #include <cstring>
  #define number s-'0'
  
  using namespace std;
  
  const int MAX_N=20;
  const int INF=0x3f3f3f3f;
  int n,m;
  int d[MAX_N][MAX_N];
  int dp[1<<MAX_N][MAX_N];
  int pre[MAX_N];
  
  void read(int &x)
  {
      char s;
      x=0;
      bool flag=0;
      while(!isdigit(s=getchar()))
          (s=='-')&&(flag=true);
      for(x=number;isdigit(s=getchar());x=x*10+number);
      (flag)&&(x=-x);
  }
  
  void write(int x)
  {
      if(x<0)
      {
          putchar('-');
          x=-x;
      }
      if(x>9) write(x/10);
      putchar(x%10+'0');
  }
  
  int min(int x, int y)
  {
      if (x<=y) return x;
      return y;
  }
  
  int main()
  {
      read(n);read(m);
      memset(d, 0x3f, sizeof(d));
      for (int i=0; i<n; i++) i[i[d]]=0;
      for (int i=0; i<m; i++)
      {
          int x,y,z;
          read(x);read(y);read(z);
          d[x][y]=z;
      }
      memset(dp, 0x3f, sizeof(dp));
      dp[(1<<n)-1][0]=0;
      for (int S=(1<<n)-2; S>=0; S--)
      {
          for (int v=0; v<n; v++)
          {
              for (int u=0; u<n; u++)
              {
                  if (!(S>>u&1))
                      dp[S][v]=min(dp[S][v], dp[S|1<<u][u]+d[v][u]);
              }
          }
      }
      write(0[0[dp]]);
      putchar('
  ');
  }

像这样针对集合的DP我们一般叫状态压缩DP

---

 Traveling by Stagecoach(POJ 2686)

• 原题如下：
  Traveling by Stagecoach

  Time Limit: 2000MS	Memory Limit: 65536K
  Total Submissions: 4482	Accepted: 1847	Special Judge

  Description

  Once upon a time, there was a traveler. 
  
  He plans to travel using stagecoaches (horse wagons). His starting point and destination are fixed, but he cannot determine his route. Your job in this problem is to write a program which determines the route for him. 
  
  There are several cities in the country, and a road network connecting them. If there is a road between two cities, one can travel by a stagecoach from one of them to the other. A coach ticket is needed for a coach ride. The number of horses is specified in each of the tickets. Of course, with more horses, the coach runs faster. 
  
  At the starting point, the traveler has a number of coach tickets. By considering these tickets and the information on the road network, you should find the best possible route that takes him to the destination in the shortest time. The usage of coach tickets should be taken into account. 
  
  The following conditions are assumed. 
  

  • A coach ride takes the traveler from one city to another directly connected by a road. In other words, on each arrival to a city, he must change the coach. 
  • Only one ticket can be used for a coach ride between two cities directly connected by a road. 
  • Each ticket can be used only once. 
  • The time needed for a coach ride is the distance between two cities divided by the number of horses. 
  • The time needed for the coach change should be ignored.

  Input

  The input consists of multiple datasets, each in the following format. The last dataset is followed by a line containing five zeros (separated by a space). 
  
  n m p a b 
  t1 t2 ... tn 
  x1 y1 z1 
  x2 y2 z2 
  ... 
  xp yp zp 
  
  Every input item in a dataset is a non-negative integer. If a line contains two or more input items, they are separated by a space. 
  
  n is the number of coach tickets. You can assume that the number of tickets is between 1 and 8. m is the number of cities in the network. You can assume that the number of cities is between 2 and 30. p is the number of roads between cities, which may be zero. 
  
  a is the city index of the starting city. b is the city index of the destination city. a is not equal to b. You can assume that all city indices in a dataset (including the above two) are between 1 and m. 
  
  The second line of a dataset gives the details of coach tickets. ti is the number of horses specified in the i-th coach ticket (1<=i<=n). You can assume that the number of horses is between 1 and 10. 
  
  The following p lines give the details of roads between cities. The i-th road connects two cities with city indices xi and yi, and has a distance zi (1<=i<=p). You can assume that the distance is between 1 and 100. 
  
  No two roads connect the same pair of cities. A road never connects a city with itself. Each road can be traveled in both directions.

  Output

  For each dataset in the input, one line should be output as specified below. An output line should not contain extra characters such as spaces. 
  
  If the traveler can reach the destination, the time needed for the best route (a route with the shortest time) should be printed. The answer should not have an error greater than 0.001. You may output any number of digits after the decimal point, provided that the above accuracy condition is satisfied. 
  
  If the traveler cannot reach the destination, the string "Impossible" should be printed. One cannot reach the destination either when there are no routes leading to the destination, or when the number of tickets is not sufficient. Note that the first letter of "Impossible" is in uppercase, while the other letters are in lowercase. 

  Sample Input

  3 4 3 1 4
  3 1 2
  1 2 10
  2 3 30
  3 4 20
  2 4 4 2 1
  3 1
  2 3 3
  1 3 3
  4 1 2
  4 2 5
  2 4 3 4 1
  5 5
  1 2 10
  2 3 10
  3 4 10
  1 2 0 1 2
  1
  8 5 10 1 5
  2 7 1 8 4 5 6 3
  1 2 5
  2 3 4
  3 4 7
  4 5 3
  1 3 25
  2 4 23
  3 5 22
  1 4 45
  2 5 51
  1 5 99
  0 0 0 0 0

  Sample Output

  30.000
  3.667
  Impossible
  Impossible
  2.856
  

  Hint

  Since the number of digits after the decimal point is not specified, the above result is not the only solution. For example, the following result is also acceptable. 
  
  

  30.0
  
  3.66667
  
  Impossible
  
  Impossible
  
  2.85595

• 题解：如果把城市看作顶点，道路看作边建图，由于有车票相关的限制，无法直接使用Dijkstra算法求解。不过，这种情况下只需要把状态作为顶点，而把状态的转移看成边建图就可以很好地避免这个问题。考虑"现在在城市v，此时还剩下的车票的集合为S"这样的状态，从这个状态出发，使用一张车票i∈S移动到相邻的城市u，就相当于转移到了"在城市u，此时还剩下的车票的集合为S{i}"这个状态。把这个转移看成一条边，那么边上的花费是(v-u间道路的长度)/t_i。按照上述的方法所构的图就可以用Dijkstra算法求解了。集合S使用状态压缩的方法表示就可以了。由于剩余的车票的集合S随着移动元素个数不断变小，因此这个图实际上一个DAG，计算DAG的最短路不需要是用Dijkstra算法，可以简单地通过DP求解。
• 代码：

  #include <cstdio>
  #include <cstdio>
  #include <cctype>
  #include <cstring>
  #include <algorithm>
  #define number s-'0'
  
  using namespace std;
  
  const int INF=0x3f3f3f3f;
  const int MAX_N=10;
  const int MAX_M=40;
  int n,m,p,a,b;
  int t[MAX_N];
  int d[MAX_M][MAX_M];
  double dp[1<<MAX_N][MAX_M];
  
  void read(int &x)
  {
      char s;
      x=0;
      bool flag=0;
      while(!isdigit(s=getchar()))
          (s=='-')&&(flag=true);
      for(x=number;isdigit(s=getchar());x=x*10+number);
      (flag)&&(x=-x);
  }
  
  double min(double x, double y)
  {
      if (x<=y) return x;
      return y;
  }
  
  int main()
  {
      read(n);read(m);read(p);read(a);read(b);
      while (n+m+p+a+b)
      {
          for (int i=0; i<n; i++) read(t[i]);
          memset(d, -1, sizeof(d));
          for (int i=0; i<p; i++)
          {
              int x, y, z;
              read(x);read(y);read(z);
              x--;y--;
              x[y[d]]=y[x[d]]=z;
          }
          for (int i=0; i<1<<n; i++)
          {
              fill(dp[i], dp[i]+m, INF);
          }
          dp[(1<<n)-1][a-1]=0;
          double res=INF;
          for (int S=(1<<n)-1; S>=0; S--)
          {
              res=min(res, dp[S][b-1]);
              for (int v=0; v<m; v++)
              {
                  for (int i=0; i<n; i++)
                   {
                       if ((S>>i)&1)
                       {
                           for (int u=0; u<m; u++)
                           {
                               if (d[v][u]>=0)
                               {
                                   dp[S&~(1<<i)][u]=min(dp[S&~(1<<i)][u],dp[S][v]+d[v][u]/(double)t[i]);
                               }
                           }
                       }
                   }
               }
          }
          if (res==INF) puts("Impossible");
          else printf("%.3f
  ", res);
          read(n);read(m);read(p);read(a);read(b);
      }
  }

---

铺砖问题

• 问题描述：给定n*m的格子，每个格子被染成黑色或者白色(.表示白色，x表示黑色)。现在要用1*2的砖块覆盖这些格子，要求块与块之间互相不重叠，且覆盖了所有白色的格子，但不覆盖任意一个黑色格子。求一共有多少种覆盖方法。输出方案数对M取余后的结果。
• 限制条件：
  1≤n≤15
  1≤m≤15
  2≤M≤10⁹
• 分析：首先考虑枚举所有的解为了不重复统计，每次从最左上方的空格处开始放置：

  #include <cstdio>
  #include <cstring>
  
  using namespace std;
  
  const int MAX_N=15;
  const int MAX_M=15;
  const int M=100007;
  int n,m;
  bool color[MAX_N][MAX_M];
  
  int rec(int i, int j, bool used[MAX_N][MAX_M]);
  
  int main(int argc, char * argv[])
  {
      scanf("%d %d
  ", &n, &m);
      for (int i=0; i<n; i++)
      {
          char c;
          for (int j=0; j<m; j++)
          {    
              scanf("%c", &c);
              if (c=='.') color[i][j]=false;
              else color[i][j]=true;
          }
          getchar();
      }
      bool used[MAX_N][MAX_M];
      memset(used, 0, sizeof(used));
      printf("%d
  ", rec(0, 0, used));
  }
  
  int rec(int i, int j, bool used[MAX_N][MAX_M])
  {
      if (j==m) return rec(i+1, 0, used);
      if (i==n) return 1;
      if (used[i][j] || color[i][j]) return rec(i, j+1, used);
      int res=0;
      used[i][j]=true;
      if (j+1<m && !used[i][j+1] && !color[i][j+1])
      {
          used[i][j+1]=true;
          res+=rec(i, j+1, used);
          used[i][j+1]=false;
      }
      if (i+1<n && !used[i+1][j] && !color[i+1][j])
      {
          used[i+1][j]=true;
          res+=rec(i, j+1, used);
          used[i+1][j]=false;
      }
      used[i][j]=false;
      return res%M;
  }

  这个方法效率不足，且递归函数的参数共有n*m*2^n*m种可能，无法使用记忆化搜索。
  自习思考之后会发现，实际上参数并没有那么多种可能，首先，由于黑色的格子不能被覆盖，因此used里对应的位置总是false，对于白色的格子，如果现在要在(i,j)位置上放置砖块，那么由于总是从最左上方的可放的格子开始放置，因此对于(i',j')<(i,j)（按字典序比较）的(i',j')总有used[i'][j']=true成立，此外，由于砖块的大小为1*2，因此对于每一列j'在满足(i',j')≥(i,j)的所有i'中，除了最小的i'之外都满足used[i'][j']=false。因此，不确定的只有每一列里还没有查询的格子中最上面的一个，共m个。从而可以把这m个格子通过状态压缩编码进行记忆化搜索，复杂度为O(n*m*2^m)。

• 代码：

  #include <cstdio>
  #include <cstring>
  
  using namespace std;
  
  const int MAX_N=15;
  const int MAX_M=15;
  const int M=100007;
  int n,m;
  bool color[MAX_N][MAX_M];
  int dp[2][1<<MAX_M];
  
  void swap(int *a[], int *b[])
  {
      int *t;
      t=*a;
      *a=*b;
      *b=t;
  }
  
  int main(int argc, char * argv[])
  {
      scanf("%d %d
  ", &n, &m);
      for (int i=0; i<n; i++)
      {
          char c;
          for (int j=0; j<m; j++)
          {    
              scanf("%c", &c);
              if (c=='.') color[i][j]=false;
              else color[i][j]=true;
          }
          getchar();
      }
      int *crt=dp[0], *next=dp[1];
      crt[0]=1;
      for (int i=n-1; i>=0; i--)
      {
          for (int j=m-1; j>=0; j--)
          {
              for (int used=0; used<1<<m; used++)
              {
                  if (used>>j&1 || color[i][j])
                  {
                      next[used]=crt[used & ~(1<<j)];
                  }
                  else
                  {
                      int res=0;
                      if (j+1<m && !(used>>(j+1)&1) && !color[i][j+1])
                      {
                          res+=crt[used | 1 << (j+1)];
                      }
                      if (i+1<n && !color[i+1][j])
                      {
                          res+=crt[used | 1 << j];
                      }
                      next[used]=res%M;
                  }
              }
              swap(&crt, &next);
          }
      }
      printf("%d
  ", crt[0]);
  }

   上面的铺砖方案，用图论的语言来说就是一个完美匹配。平面图的完美匹配的个数可以在多项式时间内求解。

---

矩阵的幂

---

斐波那契数列

• 问题描述：F₀=0 F₁=1 F_n+2=F_n+1+F_n，求第n项的值对10⁴取余后的结果
• 限制条件：0≤n≤10¹⁶
• 分析：逐项计算递推式复杂度O(n)，效率太低，如果求出通项：F_n=(1/√5){[(1+√5)/2]ⁿ-[(1-√5)/2]ⁿ}，由于式中包含无理数，无法简单求得模10⁴之后的结果。况且，一些问题很难求得通项公式。这些情况都可以不求出通项公式，而用矩阵高效地求出第n项的值。
  
  把斐波那契数列表示成矩阵就得到了下面的式子：
  
  
  
  记这个矩阵为A，则有
  
  
  
  因此只要求出Aⁿ就可以求出F_n了。Aⁿ用快速幂运算，在O(logn)时间里求出第n项的值。
• 代码：

  #include <cstdio>
  #include <cstring>
  #include <vector>
  
  using namespace std;
  
  typedef vector<int> vec;
  typedef vector<vec> mat;
  long long n;
  const int M=100007;
  
  mat mul(mat &A, mat &B)
  {
      mat C(A.size(), vec(B[0].size()));
      for (int i=0; i<A.size(); i++)
      {
          for (int k=0; k<B.size(); k++)
          {
              for (int j=0; j<B[0].size(); j++)
              {
                  C[i][j]=(C[i][j]+A[i][k]*B[k][j])%M;
              }
          }
      }
      return C;
  }
  
  mat pow(mat A, long long n)
  {
      mat B(A.size(), vec(A.size()));
      for (int i=0; i<A.size(); i++)
      {
          B[i][i]=1;
      }
      while (n>0)
      {
          if (n&1) B=mul(B, A);
          A=mul(A, A);
          n>>=1;
      }
      return B;
  }
  
  int main(int argc, char * argv[])
  {
      scanf("%d", &n);
      mat A(2, vec(2));
      A[0][0]=1;A[0][1]=1;
      A[1][0]=1;A[1][1]=0;
      A=pow(A,n);
      printf("%d
  ", A[1][0]);
  }

• 更一般地，对于m项递推式，如果记递推式为
  
  则可以把递推式写成如下矩阵形式

  
  通过计算这个矩阵的n次幂，就可以在O(m3logn)的时间内计算出第n项的值，不过，如果递推式里含有常数项则稍微复杂一些，需变成如下形式
  
  事实上，要求m项递推式的第n项的值可以不使用矩阵，而是使用初项的线性表示，用过快速幂在O(m²logn)的时间内求出答案。

---

Blocks(POJ 3734)

• 原题如下：
  Blocks

  Time Limit: 1000MS	Memory Limit: 65536K
  Total Submissions: 8004	Accepted: 3901

  Description

  Panda has received an assignment of painting a line of blocks. Since Panda is such an intelligent boy, he starts to think of a math problem of painting. Suppose there are N blocks in a line and each block can be paint red, blue, green or yellow. For some myterious reasons, Panda want both the number of red blocks and green blocks to be even numbers. Under such conditions, Panda wants to know the number of different ways to paint these blocks.

  Input

  The first line of the input contains an integer T(1≤T≤100), the number of test cases. Each of the next T lines contains an integer N(1≤N≤10^9) indicating the number of blocks.

  Output

  For each test cases, output the number of ways to paint the blocks in a single line. Since the answer may be quite large, you have to module it by 10007.

  Sample Input

  2
  1
  2

  Sample Output

  2
  6

• 题解：从最左边开始染色，设染到第i个方块为止，红绿都是偶数的方案数为a_i，红绿恰有一个是偶数的方案数为b_i，红绿都是奇数的方案数为c_i。这样，染到第i+1个方块为止，红绿都是偶数的方案数有两种可能：①到第i个方块为止红绿都是偶数，并且第i+1个方块染成了蓝色或者黄色 ② 到第i个方块为止红绿恰有一个是奇数， 并且第i+1个方块染成了奇数个对应的那种颜色。因此有递推式a_i+1=2*a_i+b_i，同样地，有b_i+1=2*a_i+2*b_i+2*c_i，c_i+1=b_i+2*c_i，把a_i，b_i，c_i的递推式用矩阵表示如下：
  
  

  因此就有：
  

• 代码：

  #include <cstdio>
  #include <cctype>
  #define number s-'0'
  #include <cstring>
  #include <vector>
  
  using namespace std;
  
  typedef vector<int> vec;
  typedef    vector<vec> mat;
  
  const int M=10007;
  int N;
  
  void read(int &x)
  {
      char s;
      x=0;
      bool flag=0;
      while (!isdigit(s=getchar()))
          (s=='-')&&(flag=true);
      for (x=number; isdigit(s=getchar());x=x*10+number);
      (flag)&&(x=-x);
  }
  
  void write(int x)
  {
      if (x<0)
      {
          putchar('-');
          x=-x;
      }
      if (x>9) write(x/10);
      putchar(x%10+'0');
  }
  
  mat mul(mat &A, mat &B)
  {
      mat C(A.size(), vec(B[0].size()));
      for (int i=0; i<A.size(); i++)
      {
          for (int j=0; j<B[0].size(); j++)
          {
              for (int k=0; k<B.size(); k++)
              {
                  C[i][j]=(C[i][j]+A[i][k]*B[k][j])%M;
              }
          }
      }
      return C;
  }
  
  mat pow(mat A, int n)
  {
      mat B(A.size(), vec(A.size()));
      for (int i=0; i<A.size(); i++) B[i][i]=1;
      while (n>0)
      {
          if (n&1) B=mul(B, A);
          A=mul(A, A);
          n>>=1;
      }
      return B;
  }
  
  int main(int argc, char * argv[])
  {
      int t;
      read(t);
      while (t>0)
      {
          t--;
          read(N);
          mat A(3, vec(3));
          A[0][0]=2;A[0][1]=1;A[0][2]=0;
          A[1][0]=2;A[1][1]=2;A[1][2]=2;
          A[2][0]=0;A[2][1]=1;A[2][2]=2;
          A=pow(A, N);
          printf("%d
  ",A[0][0]);
      }
  }

---

图中长度为k的路径的计数

• 问题描述：给定一个n个顶点边长为1的有向图的邻接矩阵，求这个图里长度为k的路径的总数。路径中同一条边允许通过多次。
• 限制条件：
  1≤n≤100
  1≤k≤10⁹
• 题解：假设从u出发，到v的长度为k的路径的总数为Gk[u][v]。G_k=G₁^k。复杂度为O(n³logn)
  

---

Matrix Power Series(POJ 3233)

• 原题如下：
  Matrix Power Series

  Time Limit: 3000MS	Memory Limit: 131072K
  Total Submissions: 28039	Accepted: 11438

  Description

  Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

  Input

  The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 10⁹) and m (m < 10⁴). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order.

  Output

  Output the elements of S modulo m in the same way as A is given.

  Sample Input

  2 2 4
  0 1
  1 1

  Sample Output

  1 2
  2 3

• 题解：构造矩阵：
  
  

  此时，令S_k=I+A+…+A^k-1，则有：
  
  通过计算这个矩阵的k次幂，就可求出A的累乘和，时间复杂度为O(n³logk)

• 代码：

  #include <cstdio>
  #include <cctype>
  #define number s-'0'
  #include <cstring>
  #include <vector>
  
  using namespace std;
  
  typedef vector<int> vec;
  typedef    vector<vec> mat;
  
  int n,k,m;
  mat A;
  
  void read(int &x)
  {
      char s;
      x=0;
      bool flag=0;
      while (!isdigit(s=getchar()))
          (s=='-')&&(flag=true);
      for (x=number; isdigit(s=getchar());x=x*10+number);
      (flag)&&(x=-x);
  }
  
  void write(int x)
  {
      if (x<0)
      {
          putchar('-');
          x=-x;
      }
      if (x>9) write(x/10);
      putchar(x%10+'0');
  }
  
  mat mul(mat &A, mat &B)
  {
      mat C(A.size(), vec(B[0].size()));
      for (int i=0; i<A.size(); i++)
      {
          for (int j=0; j<B[0].size(); j++)
          {
              for (int k=0; k<B.size(); k++)
              {
                  C[i][j]=(C[i][j]+A[i][k]*B[k][j])%m;
              }
          }
      }
      return C;
  }
  
  mat pow(mat A, int n)
  {
      mat B(A.size(), vec(A.size()));
      for (int i=0; i<A.size(); i++) B[i][i]=1;
      while (n>0)
      {
          if (n&1) B=mul(B, A);
          A=mul(A, A);
          n>>=1;
      }
      return B;
  }
  
  int main(int argc, char * argv[])
  {
      read(n);read(k);read(m);
      A=mat (n,vec(n));
      mat B(n*2, vec(n*2));
      for (int i=0; i<n; i++)
      {
          for (int j=0; j<n; j++)
          {
              read(A[i][j]);
              B[i][j]=A[i][j];
          }
          B[n+i][i]=B[n+i][n+i]=1;
      }
      B=pow(B,k+1);
      for (int i=0; i<n; i++)
      {
          for (int j=0; j<n; j++)
          {
              int a=B[n+i][j]%m;
              if (i==j) a=(a+m-1)%m;
              printf("%d%c", a, j+1==n?'
  ':' ');
          }
      }
  }

---

Minimizing maximizer(POJ 1769)

• 原题如下：
  Minimizing maximizer

  Time Limit: 5000MS		Memory Limit: 30000K
  Total Submissions: 5072		Accepted: 2057

  Description

  The company Chris Ltd. is preparing a new sorting hardware called Maximizer. Maximizer has n inputs numbered from 1 to n. Each input represents one integer. Maximizer has one output which represents the maximum value present on Maximizer's inputs. 
  
  Maximizer is implemented as a pipeline of sorters Sorter(i1, j1), ... , Sorter(ik, jk). Each sorter has n inputs and n outputs. Sorter(i, j) sorts values on inputs i, i+1,... , j in non-decreasing order and lets the other inputs pass through unchanged. The n-th output of the last sorter is the output of the Maximizer. 
  
  An intern (a former ACM contestant) observed that some sorters could be excluded from the pipeline and Maximizer would still produce the correct result. What is the length of the shortest subsequence of the given sequence of sorters in the pipeline still producing correct results for all possible combinations of input values? 
  
  Task 
  Write a program that: 
  
  reads a description of a Maximizer, i.e. the initial sequence of sorters in the pipeline, 
  computes the length of the shortest subsequence of the initial sequence of sorters still producing correct results for all possible input data, 
  writes the result. 

  Input

  The first line of the input contains two integers n and m (2 <= n <= 50000, 1 <= m <= 500000) separated by a single space. Integer n is the number of inputs and integer m is the number of sorters in the pipeline. The initial sequence of sorters is described in the next m lines. The k-th of these lines contains the parameters of the k-th sorter: two integers ik and jk (1 <= ik < jk <= n) separated by a single space.

  Output

  The output consists of only one line containing an integer equal to the length of the shortest subsequence of the initial sequence of sorters still producing correct results for all possible data.

  Sample Input

  40 6
  20 30
  1 10
  10 20
  20 30
  15 25
  30 40
  

  Sample Output

  4
  

  Hint

  Huge input data, scanf is recommended.
• 题解：首先，考虑在什么情况下可以正常工作，假设输入的第i个数是应该输出的最大值，此时在第一个满足s_k≤i≤t_k的Sorter的输出中，这个值被移动到了第t_k个位置，接下去，在一个满足s_k'≤t_k≤t_k'且k'>k的Sorter的输出中，这个值又被移动到了第t_k'个。不断重复这样的操作，如果最后可以被移动到第n个，那么就表示Maximizer可以正常工作（实际上，就是要按照线段的输入顺序，将[1, n]从左到右依次覆盖，求所需线段个数的最小值）。因此只要i=1的情况可以正常工作，那么对于任意的i都可以正常工作。不妨假设第一个数是应该输出的最大值，考虑DP：
    dp[i][j]:=到第i个Sorter为止，最大值被移动到第j个位置所需要的最短的子序列的长度（INF表示不存在这样的序列）
    dp[0][1]=0
    dp[0][j]=INF(j>1)
    dp[i+1][j]=①dp[i][j] (t_i ≠ j)
                    ②min( dp[i][j] , min{dp[i][j']|s_i≤j'≤t_i}+1 ) (t_i = j)
  由于这个DP的复杂度是O(nm)的，仍然无法在规定的时间内求出答案。但是对于(t_i ≠ j)时有dp[i+1][j]=dp[i][j]，我们可以使用同一个数组不断对自己更新：
    dp[j]:=最大值被移动到第j个位置所需要的最短的子序列的长度。(INF表示不存在这样的序列)
  初始化：
    dp[1]=0
    dp[j]=INF(j>1)
  对于每个i，这样更新：
    dp[t_i]=min(  dp[t_i], min{dp[j']|s_i≤j'≤t_i}+1 )
  这样，对于每个i都只需更新一个值就可以了，但求最小值时，最坏情况下要O(n)的时间，最后复杂度还是O(nm)。如果使用线段树来维护，就可以在O(mlogn)时间内求解了。
• 代码：

  #include <cstdio>
  #include <cctype>
  #include <algorithm>
  #define number s-'0'
  
  using namespace std;
  
  const int INF=0x3f3f3f3f;
  const int MAX_N=500005;
  const int MAX_M=500005;
  int n,m;
  int s[MAX_M], t[MAX_M];
  int dp[MAX_N+1];
  int dat[4*MAX_N];
  
  void read(int &x)
  {
      char s;
      bool flag=0;
      x=0;
      while (!isdigit(s=getchar()))
          (s=='-')&&(flag=true);
      for (x=number; isdigit(s=getchar());x=x*10+number);
      (flag)&&(x=-x);
  }
  
  void write(int x)
  {
      if (x<0)
      {
          putchar('-');
          x=-x;
      }
      if (x>9) write(x/10);
      putchar(x%10+'0');
  }
  
  void rmq_init(int k, int l, int r);
  void update(int u, int v, int k, int l, int r);
  int query(int a, int b, int k, int l, int r);
  
  int min(int x, int y)
  {
      if (x<y) return x;
      return y;
  }
  
  int main(int argc, char * argv[])
  {
      read(n);read(m);
      for (int i=0; i<m; i++)
      {
          read(s[i]);read(t[i]);
          s[i]--;t[i]--;
      }
      rmq_init(0, 0, n);
      fill(dp, dp+MAX_N, INF);
      dp[0]=0;
      update(0, 0, 0, 0, n);
      for (int i=0; i<m; i++)
      {
          int v=min(dp[t[i]], query(s[i], t[i]+1, 0, 0, n)+1);
          dp[t[i]]=v;
          update(t[i], v, 0, 0, n);
      }
      write(dp[n-1]);putchar('
  ');
  }
  
  void rmq_init(int k, int l, int r)
  {
      dat[k]=INF;
      if (r-l==1) return;
      rmq_init(k*2+1, l, (l+r)/2);
      rmq_init(k*2+2, (l+r)/2, r);
  }
  
  void update(int u, int v, int k, int l, int r)
  {
      if (r-l==1)
      {
          dat[k]=v;
          return;
      }
      else
      {
          int m=(l+r)/2;
          if (u<m) update(u, v, k*2+1, l, m);
          else update(u, v, k*2+2, m, r);
          dat[k]=min(dat[k*2+1], dat[k*2+2]);
      }
  }
  
  int query(int a, int b, int k, int l, int r)
  {
      if (r<=a || b<=l) return INF;
      if (a<=l && r<=b) return dat[k];
      else 
      {
          int m=(l+r)/2;
          int vl=query(a, b, k*2+1, l, m);
          int vr=query(a, b, k*2+2, m, r);
          return min(vl, vr);
      }
  }