矩阵乘法（三）：根据要求构造矩阵进行快速幂运算

      在应用矩阵的快速幂运算解决一些递推问题时，由于递推式不是一个直接的线性关系，这样不能直接简单地得到用于运算的矩阵，需要进行适当的构造。下面先看一道poj 上的经典题目。

【例1】matrix power series （poj 3233）

description

given a n × n matrix a and a positive integer k, find the sum s = a + a2 + a3 + … + ak.

input

the input contains exactly one test case. the first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). 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

      （1）编程思路。

      题目的意思是：输入n*n矩阵a，求s=a + a^2 + a^3 + … + a^k的结果（两个矩阵相加就是对应位置分别相加），输出的数据mod m。

      显然不能将各个矩阵的2~k次幂均计算出来，会超时的。

      可以构造出一个2n*2n的矩阵p，采用分块的方法表示如下，，其中 a 为原矩阵，i 为单位矩阵，o 为0矩阵。

    

      

        我们发现 p^k+1 右上角那一个分块n*n矩阵正是要求的 a+a²+...+a^k 。

　　于是我们构造出 p矩阵，然后对它求矩阵快速幂，最后减去一个单位阵即可得到结果。

      （2）源程序。

#include <stdio.h>
#include <string.h>
struct matrix
{
      int mat[61][61]; // 存储矩阵中各元素
};
matrix matmul(matrix a ,matrix b,int n,int m)
{
      matrix c;
      memset(c.mat,0,sizeof(c.mat));
      int i,j,k;
      for (i = 1; i<=n ; i++)
          for (j=1 ;j<=n ; j++)
              for (k = 1 ;k<=n ;k++)
              {
                    c.mat[i][j]=(c.mat[i][j]+a.mat[i][k] * b.mat[k][j]) % m;
              }
      return c;
}
matrix quickmatpow(matrix a ,int n,int b,int m) // n阶矩阵a快速b次幂
{
      matrix c;
      memset(c.mat ,0 ,sizeof(c.mat));
      int i;
      for (i = 1 ;i <= n ;i++)
      c.mat[i][i] = 1;
      while (b!=0)
      {
            if (b & 1)
                 c = matmul(c ,a ,n,m); // c=c*a;
            a = matmul(a ,a ,n,m); // a=a*a
            b /= 2;
      }
      return c;
}
int main()
{
      int n,k,m,i,j;
      matrix p ;
      scanf("%d%d%d",&n,&k,&m);
      memset(p.mat,0,sizeof(p.mat));
      for (i=1;i<=n;i++)
          for (j=1;j<=n;j++)
              scanf("%d",&p.mat[i][j]);
      for (i=1;i<=n;i++)
      {
               p.mat[i][n+i]=1;        // 右上角的单位矩阵
               p.mat[n+i][n+i]=1;    // 右下角的单位矩阵
      }
      p = quickmatpow(p,2*n,k+1,m);
      for (i=1;i<=n;i++)
      p.mat[i][i+n]--;               // 减去单位矩阵
      for (i=1;i<=n;i++)
      {
            for (j=n+1;j<=2*n;j++)
                 printf("%d ",(p.mat[i][j]+m)%m);
            printf("\n");
      }
      return 0;
}

      由于构造的p矩阵中含有大量的0元素，因此可以考虑在矩阵相乘时进行优化，优化的方法是如果是0元素就不进行对应元素运算。优化后的矩阵相乘函数如下：

matrix matmul(matrix a ,matrix b,int n,int m)
{
      matrix c;
      memset(c.mat,0,sizeof(c.mat));
      int i,j,k;
      for (k = 1; k<=n ; k++)
          for (i=1 ;i<=n ; i++)
             if (a.mat[i][k]!=0)
                 for (j = 1 ;j<=n ;j++)
                     c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j]) % m;
      return c;
}

【例2】另类斐波那契数列。

      已知斐波那契数列为：f(0) = 1，f(1) = 1，f(n) = f(n - 1) + f(n - 2)  (n >= 2)。现定义一个新的类斐波那契数列： a(0) = 1，a(1) = 1，a(n) = x * a(n - 1) + y * a(n - 2) (n >= 2)。

      输入x、y和n，求新定义数列的前n项的平方和，即s(n)=a(0)^2+a(1)^2+....+a(n)^2。

      （1）编程思路。

      根据数列递推式不能直观地得到用于递推的矩阵，需要进行分析、构造。 

      因为  s[ n ]  = s[ n -1 ]  + a[n]^2

               a[n]=x*a[n-1]+y*a[n-2]

      所以  s[n]=s[n-1]+x^2*a[n-1]^2+y^2*a[n-2]^2+2xya[n-1]*a[n-2]

      等式的右边有四项，有变化的是s[ n-1]，a[n-1]^2，a[n-2]^2，a[n-1]*a[n-2]，各自的系数1，x^2，y^2，2xy一经输入就确定了，不会再变化。

      而有变化的四项s[ n-1]，a[n-1]^2，a[n-2]^2，a[n-1]*a[n-2]向前推进一步应该为

s[ n]，a[n]^2，a[n-1]^2，a[n]*a[n-1]。

      这样，可以构造一个4*4的矩阵p。

          

 

     为什么这样构造呢？是因为   （仔细体会一下哟!）

      

 

      构造好矩阵p后，就可以采用矩阵快速幂运算求s(n)了。

       

 

      （2）源程序。

#include <stdio.h>
#include <string.h>
#define mod 10007
struct matrix
{
      int mat[5][5]; // 存储矩阵中各元素
};
matrix matmul(matrix a ,matrix b,int n)
{
      matrix c;
      memset(c.mat,0,sizeof(c.mat));
      int i,j,k;
      for (k = 1; k<=n ; k++)
          for (i=1 ;i<=n ; i++)
              if (a.mat[i][k]!=0)
                  for (j = 1 ;j<=n ;j++)
                      c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j]) % mod;
      return c;
}
matrix quickmatpow(matrix a ,int n,int b) // n阶矩阵a快速b次幂
{
      matrix c;
      memset(c.mat ,0 ,sizeof(c.mat));
      int i;
      for (i = 1 ;i <= n ;i++)
             c.mat[i][i] = 1;
      while (b!=0)
      {
            if (b & 1)
                c = matmul(c ,a ,n); // c=c*a;
            a = matmul(a ,a ,n); // a=a*a
            b /= 2;
      }
      return c;
}
int main()
{
       int n,x,y,ans;
       matrix p;
       while(scanf("%d%d%d" ,&n ,&x ,&y)!=eof)
       {
            x = x%mod;
            y = y%mod ;
            if (n==2)
                 printf("%d\n" ,(x*x%mod+y*y%mod+2*x*y%mod+2)%mod) ;
            else
            {
                 memset(p.mat,0,sizeof(p.mat));
                 p.mat[1][1]=p.mat[3][2]=1;
                 p.mat[1][2]=p.mat[2][2]=(x*x)%mod;
                 p.mat[1][3]=p.mat[2][3]=(y*y)%mod;
                 p.mat[1][4]=p.mat[2][4]=(2*x*y)%mod;
                 p.mat[4][2]=x;
                 p.mat[4][4]=y;
                 p = quickmatpow(p,4,n-1);
                 ans=(p.mat[1][1]*2%mod+p.mat[1][2]+p.mat[1][3]+p.mat[1][4])%mod;
                 printf("%d\n" ,ans);
            }
       }
       return 0;
}

       将此源程序提交给 hdu 3306 “another kind of fibonacci”，可以accepted。   

【例3】又一个非线性递推数列。

      设有数列f的定义如下：f[1] =a，f[2]=b，  f(n)=f(n-1)+f(n-2)*2+n^4  （n>=3）。

      输入a、b和n的值， 求f[n] mod 2147493647的结果。

      （1）编程思路。

      通过本题可以更进一步思考用于递推的矩阵的构造方法。

      因为  f(n)=f(n-1)+f(n-2)*2+n^4

               f(n-1)=f(n-2)+f(n-3)*2+(n-1)^4

      因此，构造矩阵时的关键是找到怎样由(n-1)^4到n^4的递推。

      又由于 n^4=(n-1+1)^4=(n-1)^4 + 4 (n-1)^3 + 6 (n-1)^2 + 4 (n-1)^1 +1，因此，需要通过构造

矩阵维护好f(n-1)、f(n-2)、(n-1)^4、(n-1)^3、(n-1)^2、(n-1)^1和(n-1)^0这7个值。

     可以构造p矩阵如下：

        

 

       为什么这样构造呢？是因为 （好好体会哟！）

 

      构造好矩阵p后，就可以采用矩阵快速幂运算求f(n)了。

     

 

      （2）源程序。

#include <stdio.h>
#include <string.h>
#define mod 2147493647
struct matrix
{
      __int64 mat[8][8]; // 存储矩阵中各元素
};
matrix matmul(matrix a ,matrix b,int n)
{
      matrix c;
      memset(c.mat,0,sizeof(c.mat));
      int i,j,k;
      for (k = 1; k<=n ; k++)
          for (i=1 ;i<=n ; i++)
              if (a.mat[i][k]!=0)
                  for (j = 1 ;j<=n ;j++)
                      c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j]) % mod;
      return c;
}
matrix quickmatpow(matrix a ,int n,int b) // n阶矩阵a快速b次幂
{
      matrix c;
      memset(c.mat ,0 ,sizeof(c.mat));
      int i;
      for (i = 1 ;i <= n ;i++)
           c.mat[i][i] = 1;
      while (b!=0)
      {
           if (b & 1)
               c = matmul(c ,a ,n); // c=c*a;
           a = matmul(a ,a ,n); // a=a*a
           b /= 2;
      }
      return c;
}
int main()
{
      int t,n,a,b;
      __int64 ans;
      matrix p;
      scanf("%d" ,&t);
      while(t--)
      {
           scanf("%d%d%d" ,&n ,&a ,&b);
           if (n==1)
               printf("%d\n",a);
           else if (n==2)
               printf("%d\n",b);
           else
           {
                  memset(p.mat,0,sizeof(p.mat));
                  p.mat[1][1]=p.mat[1][3]=p.mat[1][7]=1;
                  p.mat[2][1]=p.mat[3][3]=p.mat[3][7]=1;
                  p.mat[4][4]=p.mat[4][7]=p.mat[5][5]=1;
                  p.mat[5][7]=p.mat[6][6]=p.mat[6][7]=p.mat[7][7]=1;
                  p.mat[1][2]=p.mat[5][6]=2;
                  p.mat[4][5]=p.mat[4][6]=3;
                  p.mat[1][4]=p.mat[1][6]=p.mat[3][4]=p.mat[3][6]=4;
                  p.mat[1][5]=p.mat[3][5]=6;
                  p = quickmatpow(p,7,n-2);
                  ans=(b*p.mat[1][1])% mod;
                  ans=(ans+a*p.mat[1][2]) % mod;
                  ans=(ans+16*p.mat[1][3]) % mod;
                  ans=(ans+8*p.mat[1][4]) % mod;
                  ans=(ans+4*p.mat[1][5]) % mod;
                  ans=(ans+2*p.mat[1][6]) % mod;
                  ans=(ans+p.mat[1][7])% mod;
                  printf("%i64d\n" ,ans);
           }
      }
      return 0;
}

      将此源程序提交给 hdu 5950 “recursive sequence”，可以accepted。