本贴声明:
关于这道题的基本解法, 我在之前曾经发表过, 以动态规划的方式在O(N)的时间复杂度内求解, 但对于数据规模为10^9的数据而已, O(N)显然是不够的, 当时我受困良久. 但幸运的是, 某网友给了我一个万分有用的建议, 以矩阵的方式的进行求解. 当我实现以后, 我发现这是一个O( log(2)(N) )[注: 以2为底N的对数 ]的算法, 速度之快, 毋庸置疑.
这道题的矩阵解法, 其核心理论依然依靠于基本的动态规划解法, 然后再配合上线性代数中的矩阵加以优化, 如果读者不了解此题的动态解法, 那么以下解释无疑相当于看天书, 所以, 我建议读者在阅读下文之前, 先阅读我以前写的叠骰子 --- 动态规划解法的解题报告,然后, 如果你没有线性代数的基础知识, 也是不行的,最起码的, 你应该对矩阵的乘法规则有一定程度的了解, 最后, 你还应该熟悉快速幂的原理, 如果你具备上述的三项基础知识, 那么我们可以接着看下去了.
线性组合的奇妙:
在讲如何用矩阵来优化动态规划算法之前, 我先引入一个新的例子来表达矩阵乘法的奇妙之处. 相信大家都曾经接触过斐波那契数列,数列定义如下:
1) F(0) = 0, F(1) = 1;
2) F(n) = F(n-1) + F(n-2); [ n > 1 ]
相信但凡学过递归的同学都接触过该数列, 但是我们一般的解法都是以带记忆数组的递归来求解, 不过, 这个很特别的数列确是可以以矩阵来求解的, 看看下面的矩阵递归式:
式子十分神奇吧!! 通过对矩阵乘法的巧妙应用, 我们可以写成如此精致的递归式 . 仔细分析, 其实中间的矩阵的列向量是十分特别的, 它以近乎完美的方式对左边矩阵中的点值进行了有机性的取舍, 从而演变出了右边的矩阵. 而右边的矩阵与左边的矩阵又具有相似的结构, 所以, 递归式可以无限地演变下去. 上述的式子还可以继续精简 :
通过以上演变, 我们把斐波那契数列的求解变成了求矩阵的幂, 这时候, 矩阵快速幂的算法就派上用场了, 最后, 我们便可以在O( log(2)(N) )的复杂度内求解数列.
那么, 我们如何在叠骰子里应用上述算法 ? 在本题的动态规划解法中, 我曾经说过, 以dp[ i ][ j ]来表示在高度为i的情况下, 骰子柱顶面点数为j的总方案数, 然后, dp[ i ][ j ]的值将等于前一高度所有方案的选择性累加, 一说到选择性累加, 是不是十分类似与上述例子中的那个只包含0和1的矩阵 ? 在上述的斐波那契例子中, 那个只包含0和1的矩阵就是这样通过列向量的0或1对左边矩阵进行选择累加, 从而得出右边矩阵的. 那么, 我们也可以通过类似的结构来完成对骰子方案的选择性累加.
接下来, 我们需要沿用在DP解法中提到的冲突数组, 因为冲突数组实际上就是我们的选择性累加, 在这里, 我们将它变为冲突矩阵, 并且更正一下对它的定义 :
1 . Complict[ i ][ j ] 表示: 有一个点数为i的面, 它的反面与点数为j的面存在冲突.
上述定义意味着, 如果1和2存在冲突, 那么Compact[4][2] 和 Compact[5][1] 都为0 , 因为4的反面是1, 5的反面是2.
然后, 我们一个单行矩阵Count[1][ j ]来记录高度为N时, 顶面为j的总方案数.
好了, 现在便是整个算法的核心部分, 矩阵Count记录了当前某高度下的各个面的总方案数, 而矩阵Complict的列向量描述了"选择", 这已经跟斐波那契例子中的结构相当类似了, 从而我们可以证明:
Count * Complict = newCount;
如果 Count 记录了高度为N的各个面的总方案数;
那么 newCount 记录了高度为N+1的各个面的总方案数;
具体的证明我就不写出来了, 大家可以类比一下斐波那契列的例子....... 啊!!!多么完美的一个递归式啊!!好 , 现在我们可以来考虑一下矩阵Count的初始化值, 不用多说, 但N=1时, 矩阵Count中的每一个元素都等于1, 这一点, 与我在DP解法中的初始化时一样的 !
所以, 当高度为N时, 式子是这样的 :
Count = Count * Complict * Complict *........*Complict ; ( Complict被乘了N-1次)
那么就相当于求Complict的N-1次幂了.
Count = Count *( Complict^( N-1 ) );
最后, 我们把Count里面的所有方案数累加, 就得到了总方案数, 别忘啦, 还要乘以4的N次幂 !!
对于4的N次幂的求法, 我们依然采用快速幂的算法, 所以, 整个算法就成功地在O( log(2)(N) )内求解了, 对于10^9的数据, 连1秒都不用就出结果 !!! 不说废话, 上代码 !
#include <iostream>
#include <vector>
using namespace std;
typedef unsigned int uint;
typedef uint SizeType;// ...表示大小
typedef uint IndexType;// ...表示下标
typedef long long ValueType;// ...表示值
typedef vector< ValueType > RowType;// ...表示矩阵的一行
typedef vector< RowType > ContainerType;// ...矩阵容器
#define MOD 1000000007
// ...矩阵类的实现
class Matrix{
public:
Matrix( SizeType Row, SizeType Col, ValueType Init_Val=0 ):m_row(Row),m_col(Col){
m_mtx = new ContainerType( m_row, RowType(m_col,Init_Val) );
}
Matrix( const Matrix& Copyright ):m_row(Copyright.m_row),m_col(Copyright.m_col){
m_mtx = new ContainerType( *Copyright.m_mtx );
}
inline SizeType GetRow()const{ return m_row; }
inline SizeType GetCol()const{ return m_col; }
// ...矩阵赋值
const Matrix& operator=( const Matrix& right ){
*m_mtx = *right.m_mtx;
m_row = right.m_row;
m_col = right.m_col;
}
// ...获取矩阵中第rI行cI列的元素
inline const ValueType& operator()( IndexType rI,IndexType cI )const{
return (*m_mtx)[rI][cI];
}
inline ValueType& operator()( IndexType rI, IndexType cI){
return const_cast<ValueType&>( static_cast<const Matrix&>(*this)(rI,cI) );
}
// ...获得一个n阶的单位矩阵
static Matrix GetIdentityMatrix( SizeType nDimension )
{
Matrix I( nDimension,nDimension,0);
for( IndexType i = 0; i < nDimension; ++i)
I(i,i) = 1;
return I;
}
// ...矩阵乘法
const Matrix operator*( const Matrix& right )const{
Matrix Res( m_row, right.m_col,0);
if( m_col == right.m_row ){
const Matrix& left = *this;
for( IndexType i = 0; i < m_row; ++i)
for( IndexType j = 0; j < right.m_col; ++j)
for( IndexType t = 0; t < m_col; ++t )
Res(i,j) = (Res(i,j)+left(i,t)*right(t,j)%MOD)%MOD;
}
return Res;
}
const Matrix& operator*=( const Matrix& right ){return *this = (*this)*right;}
// ...矩阵快速幂
const Matrix operator^( uint N ) const{
Matrix Res = GetIdentityMatrix( m_row ),t = *this;
while( N ){
if( N&1 ) Res*=t;
t*=t;
N>>=1;
}
return Res;
}
const Matrix& operator^=(uint N){ return *this = (*this)^N; }
protected:
SizeType m_row, m_col;// ...行数, 列数
ContainerType *m_mtx;// ...容器
};
// ...常数快速幂
ValueType Pow( ValueType a, uint N ) {
ValueType Res = 1;
while( N ){
if( N&1 ) Res=Res*a%MOD;
a=a*a%MOD;
N>>=1;
}
return Res;
}
int main(int argc, char** argv) {
Matrix Complict(6,6,1); // ...冲突组合记录
Matrix Count( 1,6,1 );// ...总数记录
const IndexType Parner[6] = {3,4,5,0,1,2}; // ...反面记录
uint N, M, s1, s2;
cin >> N >> M;
for( uint i = 0; i < M; ++i){
cin >> s1 >> s2;
Complict( Parner[s1],s2 ) = Complict( Parner[s2],s1 ) = 0;
}
Complict^=(N-1);
Count*=Complict;
ValueType sum = 0;
for( IndexType i = 0; i < Count.GetCol(); ++i)
sum = (sum+Count(0,i))%MOD;
cout << sum*Pow(4,N)%MOD;
return 0;
}