大佬博客：https://blog.csdn.net/baidu_23081367/article/details/52347256

代码：
 

const int mat_size = 5;//矩阵大小，需要乘以2，为了&运算的时候需要二倍的矩阵大小

const int mod=1e9+7;

struct Matrix
{
	long long a[mat_size][mat_size];
	int x, y;//长宽
	Matrix()	//返回0矩阵
	{
		memset(a,0,sizeof(a));
	}
	Matrix(int x,int y)//返回0矩阵，并且x,y赋值
	{
		this->x = x;
		this->y = y;
		memset(a, 0,sizeof(a));
	}
 
	Matrix(int n)	//返回n*n的【单位矩阵】
	{
		this->x=n;
		this->y=n;
		memset(a,0,sizeof(a));
		for (int i = 0; i <n;++i)	a[i][i]=1;
	}
	Matrix operator * (const Matrix &B)//矩阵乘法
	{
		Matrix tmp;
		for (int i = 0; i < x; ++ i)
			for (int j = 0; j < B.y; ++ j)
			{
				tmp.a[i][j] = 0;
				for (int k = 0; k < y; ++ k)
				{
					tmp.a[i][j] = (tmp.a[i][j] + a[i][k] * B.a[k][j] % mod) % mod;
				}
			}
		tmp.x = x;
		tmp.y=B.y;
		return tmp;
	}
	Matrix operator ^ (int b)//矩阵A的b次方
	{
		Matrix ret = Matrix(x);  
		Matrix A = *this;
		while( b )  
		{  
			if( b & 1 )	ret = ret * A ;  
			b >>= 1 ;  
			A = A * A ;  
		}  
		return ret ;  
	}
	Matrix operator & (int b)//A^0 + A^1+A^2+A^3+++A^n，其中A是矩阵。最后返回的就是一个矩阵
	{
		Matrix ret = *this;
		for (int i = ret.x; i < ret.x * 2; ++ i)	
		{
			ret.a[i-ret.x][i]= 1;
			ret.a[i][i] = 1;
		}
		ret.x <<= 1;
		ret.y <<= 1;
		//pg(ret);
		ret = ret^b;
		ret.x >>= 1;
		ret.y >>= 1;
		for (int i = 0; i < ret.x; ++ i)	
			for (int j = 0; j < ret.y; ++ j)
				(ret.a[i][j] += ret.a[i][j + ret.x])%=mod;
		return ret;
	}
	void pg(Matrix A)
	{
		for (int i = 0; i <A.x; ++i)
		{
			for (int j = 0; j < A.y;++j)	cout<<A.a[i][j]<<" ";cout<<endl;
		}
		cout<<endl;
 
	}
};