HDU2294_DP+矩阵加速(实在妙题)

/*
*State: 2294    609MS    944K    4236 B    C++
*题目大意:
*        有k种珍珠,每种有n个,然后要求组合成长度为1~n的项链的总数。
*        (项链的长度为珍珠的个数),并要求项链中至少含有k种珍珠。
*解题思路:
*        复杂的组合题,至今不理解其递推式,但是知道dp的状态转移表达式为
*        f[i][j] = f[i - 1][j - 1] * (k - j + 1) + f[i - 1][j] * j;
*        f[i][j]下标的意思是有i颗珍珠时,含有j种能组合的个数。然后有了
*        这个状态转移,可以得出结果为:f[1][k] + f[2][k] + …… + f[n][k],
*        然后发现这个状态转移够囧,n<1 000 000 000,那么一般根据二维dp要
*        两维来实现,结果复杂度就变成了O(nk),超级计算机也难以接受。
*        考虑用别的方式来实现dp,猥琐地搜了解题报告后发现,原来可以用矩阵
*        加速,有一步跳跃比较高,就是一个矩阵里面把其中一维全部包括了。
*        F[i] = F[i - 1] * A;其中A为转移矩阵。然后F[i]这个矩阵要设定为:
*        f[i-1][0] f[i-1][1] f[i-1][2] …… f[i-1][k]
*        0         0         0            0
*       0         0         0            0
*        …………
*        0         0         0            0
*        这个矩阵有k+1维,为什么是k+1维,自己想。
*        转移矩阵A为
*        0 k 0   0 …… 0 0 0 
*        0 1 k-1 0 …… 0 0 0 
*        0 0 2   0 …… 0 0 0
*        ……
*        0 0 0   0 …… 0 0 k
*        为什么转移矩阵是如此的?把上面的F[i]跟转移矩阵相乘可发现跟状态转移
*        方程一致。
*        
*        有了矩阵的表达式,那么要求结果必须先求F[1] + F[2] + F[3] + …… + F[k]
*        转化为F[1] + F[1]*A + F[1]*A^2 + …… + F[1]*A^(k-1).提取F[1],之后就是
*        二分求sum(A^n)了。
*题目困惑:
*        题意看了好久,感觉题目表示得不好,表示一条项链必须含有k钟珍珠的时候。
*        to show his wealth, every kind of pendant must be made of K pearls.
*        好像是在说需要k个珍珠一样,不是应该加个 K kind of pearls么?
*        最后理清思路了,还是wa了3次,因为最后的模溢出了,还是自己走投无路给自己
*        乱出数据得出来的,这个在实在没办法的时候,不失为一种方法。
*/
View Code
  1 #include <stdio.h>
  2 #include <string.h>
  3 #include <iostream>
  4 
  5 #define MAX_DIMENSION 35 
  6 using namespace std;
  7 
  8 typedef __int64 MATRIX_TYPE;
  9 typedef __int64 MAX_INT_TYPE; 
 10 typedef MATRIX_TYPE Matrix[MAX_DIMENSION][MAX_DIMENSION];
 11 int ndim=MAX_DIMENSION;
 12 int mod=1234567891;//mod
 13 
 14 void m_zero(Matrix  x)
 15 {
 16     memset(x, 0, sizeof(MATRIX_TYPE)*MAX_DIMENSION*MAX_DIMENSION);
 17 }
 18 
 19 void m_one(Matrix  x)
 20 {
 21     int i;
 22     m_zero(x);
 23     for(i=0;i<ndim;++i)x[i][i]=1;
 24 }
 25 
 26 void m_copy(Matrix  src,Matrix  dest)
 27 {
 28     memcpy(dest,src, sizeof(MATRIX_TYPE)*MAX_DIMENSION*MAX_DIMENSION);
 29 }
 30 
 31 //z=x+y;
 32 void m_add(Matrix  x,Matrix  y,Matrix  z)
 33 {
 34     int i,j;
 35     for(i=0;i<ndim;i++)
 36         for(j=0;j<ndim;j++)
 37             if(mod<=1)z[i][j]=x[i][j]+y[i][j];
 38             else z[i][j]=(x[i][j]+(MAX_INT_TYPE)y[i][j])%mod;//module
 39 }
 40 
 41 
 42 //c=a*b
 43 void m_multiple(Matrix  a,Matrix b,Matrix c)
 44 {
 45     int i,j,k;
 46     MAX_INT_TYPE t;
 47 
 48     for(i=0;i<ndim;i++)
 49         for(j=0;j<ndim;j++)
 50         {
 51             t=0;
 52             if(mod<=1)
 53                 for(k=0;k<ndim;k++) t+=a[i][k]*b[k][j];//module
 54             else
 55                 for(k=0;k<ndim;k++){
 56                     t+=(a[i][k]*(MAX_INT_TYPE)b[k][j])%mod;
 57                     t%=mod;
 58                 }//module
 59             c[i][j]=t;
 60         }
 61 }
 62 
 63 void m_pow_with_saved(Matrix  x_p[],unsigned int n, Matrix y)
 64 {
 65     Matrix temp;
 66     m_one(y);
 67     for(int k=1;n;++k,n>>=1)
 68     {
 69         if ((n & 1) != 0)
 70         {
 71             m_multiple(x_p[k],y,temp);
 72             m_copy(temp,y);
 73         }
 74     }
 75 }
 76 
 77 void m_pow_sum1(Matrix  x_p[],unsigned int n, Matrix y)
 78 {
 79     m_zero(y);
 80     if(n==0)return;
 81     if(n==1) m_copy(x_p[1],y);
 82     else
 83     {
 84         Matrix temp;
 85         //calculate x^1+...+^(n/2)
 86         m_pow_sum1(x_p,n>>1,temp);
 87         //add to y
 88         m_add(temp,y,y);
 89         //calculate x^(1+n/2)+...+x^n
 90         Matrix temp2;
 91         m_pow_with_saved(x_p,n>>1,temp2);
 92         Matrix temp3;
 93         m_multiple(temp,temp2,temp3);
 94         //add to y
 95         m_add(temp3,y,y);
 96         if(n&1)
 97         {
 98             //calculate x^(n-1)
 99             m_multiple(temp2,temp2,temp3);
100             //calculate x^n
101             m_multiple(temp3,x_p[1],temp2);
102             //add x^n
103             m_add(temp2,y,y);
104         }
105     }
106 
107 }
108 
109 void m_pow_sum(Matrix x, unsigned int n, Matrix y)
110 {
111     //calculate x^1 x^2 x^4 ... x^logn
112     unsigned int i=0,logn,n2=n;
113     for(logn=0,n2=n;n2;logn++,n2 >>= 1);
114     Matrix *pow_arry=new Matrix[logn==0?2:(logn+1)];
115     m_one(pow_arry[0]);
116     m_copy(x,pow_arry[1]);
117     for(i=1;i<logn;i++)
118     {
119         m_multiple(pow_arry[i],pow_arry[i],pow_arry[i+1]);
120     }
121 
122     m_pow_sum1(pow_arry,n,y);
123     delete []pow_arry;
124 }
125 
126 void view_mat(Matrix a, int n)
127 {
128     for(int i = 0; i < n; i++)
129     {
130         for(int j = 0; j < n; j++)
131             cout << a[i][j] << " ";
132         cout << endl;
133     }
134 }
135 
136 int main(void)
137 {
138 #ifndef ONLINE_JUDGE
139     freopen("in.txt", "r", stdin);
140 #endif
141 
142     int cas;
143     scanf("%d", &cas);
144     while(cas--)
145     {
146         int n, k;
147         scanf("%d %d", &n, &k);
148         Matrix a;
149         ndim = k + 1;
150         m_zero(a);
151         for(int i = 0; i <= k; i++)
152         {
153             if(i == 0)
154                 continue;
155             a[i][i] = i;
156             a[i - 1][i] = k - i + 1;
157         }    
158         Matrix asum, res, addsum, f1, ans;
159         m_one(res);
160         m_pow_sum(a, n - 1, asum);
161         m_add(res, asum, addsum);
162         m_zero(f1);
163         f1[0][1] = k;
164         m_multiple(f1, addsum, ans);
165         printf("%d\n", ((ans[0][k] % mod) + (__int64)mod) % mod);
166     }
167     return 0;
168 }

这道题目把dp的弱性跟矩阵的优美强悍地给表现出来了。

原文地址:https://www.cnblogs.com/cchun/p/2619508.html