bzoj 4162 shlw loves matrix II - 行列式 - 矩阵乘法 - 高斯消元

题目大意：
给一个n * n的矩阵A，求其k次方。$n\le 50,k\le 2^{10000}$。
题解：
考虑将$A^k$视为某个数列的第k项$h_k$（严格意义上是矩阵列？）
尝试改造$h$的递推式。
考虑：$$F(x)=|A-xI|$$即A减去若干倍的单位矩阵求行列式。
以3*3的矩阵为例，它张这个样子：
$$F(x)=\left|\begin{array}{ccc}{A}_{1,1}-x& A_{1,2}& A_{1,3}\\ A_{2,1}& A_{2,2}-x& A_{2,3}\\ A_{3,1}& A_{3,2}& A_{3,3}-x\end{array}\right|$$
显然其会是一个关于x的n次多项式。
那么根据定义有：$F(A)=0$（啥你问我为啥带入的不是实数而是个矩阵？其实严格的说法是，把F的系数copy一遍弄出一个新的以矩阵为变量的函数F’，然后根据某个C开头的定理F’(A)=0，不过直观"理解"$F(A)=|A-IA|=0$）
总之把这个多项式的系数插值出来，可以知道：
$$\sum _{i=0}^n{f}_i{A}^i=0$$
那么我要求第k项：
$$\sum _{i=0}^n{f}_i{A}^{i+k-n}=\sum _{i=0}^n{f}_{n-i}{A}^{k-i}=0$$
或者说：
$$\sum _{i=0}^n{f}_{n-i}{h}_{k-i}=0$$
因此可以通过$h_{k-n},\dots ,h_{k-1}$算出$h_k$，套用常系数线性递推即可。
复杂度$O\left(n^4+n^{2}lgk\right)$，复杂度瓶颈在于算n次行列式和预处理$h_0,...,h_{n-1}$。

#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define Rep(i,v) rep(i,0,(int)v.size()-1)
#define lint long long
#define mod 1000000007
#define ull unsigned lint
#define db long double
#define pb push_back
#define mp make_pair
#define fir first
#define sec second
#define gc getchar()
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
using namespace std;
typedef pair<int,int> pii;
typedef set<int>::iterator sit;
inline int inn()
{
	int x,ch;while((ch=gc)<'0'||ch>'9');
	x=ch^'0';while((ch=gc)>='0'&&ch<='9')
		x=(x<<1)+(x<<3)+(ch^'0');return x;
}
const int N=53;int v[N][N],vs[N][N],a[N][N],b[N];
inline int fast_pow(int x,int k,int ans=1) { for(;k;k>>=1,x=(lint)x*x%mod) (k&1)?ans=(lint)ans*x%mod:0;return ans; }
inline int calc(int n,int x)
{
	rep(i,1,n) memcpy(vs[i],v[i],sizeof(int)*(n+1));
	rep(i,1,n) v[i][i]-=x,(v[i][i]<0?v[i][i]+=mod:0);
	int det=1;
	rep(i,1,n)
	{
		int x=i;rep(j,i,n) if(v[j][i]) x=j;
		if(i^x) swap(v[i],v[x]),det=(det?mod-det:0);
		if(!v[i][i]) { det=0;break; }x=fast_pow(v[i][i],mod-2);
		rep(j,i,n) v[i][j]=(lint)v[i][j]*x%mod;det=(lint)det*x%mod;
		rep(j,i+1,n) { x=v[j][i];rep(k,i,n) v[j][k]-=(lint)x*v[i][k]%mod,(v[j][k]<0?v[j][k]+=mod:0); }
	}
	rep(i,1,n) memcpy(v[i],vs[i],sizeof(int)*(n+1));
	return det?fast_pow(det,mod-2):0;
}
inline int gauss(int (*a)[N],int *b,int n)
{
	rep(i,0,n)
	{
		int x=i;rep(j,i,n) if(a[j][i]) x=j;
		swap(a[i],a[x]),swap(b[i],b[x]),x=fast_pow(a[i][i],mod-2);
		rep(j,i,n) a[i][j]=(lint)a[i][j]*x%mod;b[i]=(lint)b[i]*x%mod;
		rep(j,0,n) if(i^j)
		{
			x=a[j][i];if(!x) continue;
			rep(k,i,n) a[j][k]-=(lint)x*a[i][k]%mod,(a[j][k]<0?a[j][k]+=mod:0);
			b[j]-=(lint)x*b[i]%mod,(b[j]<0?b[j]+=mod:0);
		}
	}
	return 0;
}
struct matrix{
#define P(k) (lint)a.v[i][k]*b.v[k][j]
	int v[N][N],n;matrix(int _n=0) { init(_n); }
	inline int init(int _n=0) { n=_n;rep(i,1,n) memset(v[i],0,sizeof(int)*(n+1));return 0; }
	inline matrix operator*(const matrix &b)const
	{
		const matrix &a=*this;matrix c(n);
		rep(i,1,n) rep(k,1,n) rep(j,1,n) c.v[i][j]=(c.v[i][j]+P(k))%mod;
		return c;
	}
	inline matrix& operator*=(const matrix &b) { return (*this)=(*this)*b; }
	inline matrix operator+(const matrix &b)const
	{
		const matrix &a=*this;matrix c(n);
		rep(i,1,n) rep(j,1,n) c.v[i][j]=a.v[i][j]+b.v[i][j],(c.v[i][j]>=mod?c.v[i][j]-=mod:0);
		return c;
	}
	inline matrix& operator+=(const matrix &b) { return (*this)=(*this)+b; }
	inline matrix operator*(int x)const
	{
		const matrix &a=*this;matrix c(n);
		rep(i,1,n) rep(j,1,n) c.v[i][j]=(lint)a.v[i][j]*x%mod;
		return c;
	}
	inline matrix& operator*=(int x) { return (*this)=(*this)*x; }
	inline int print()const { rep(i,1,n) rep(j,1,n) printf("%d%c",v[i][j]," \n"[j==n]);return 0; }
}A,B,ans;
struct poly{
	int v[N<<1],n;
	poly() { memset(v,0,sizeof v),n=0; }
	inline poly operator*(const poly &b)const
	{
		const poly &a=*this;poly c;
		c.n=a.n+b.n;
		rep(i,0,a.n) rep(j,0,b.n)
			c.v[i+j]=(c.v[i+j]+(lint)a.v[i]*b.v[j])%mod;
		return c;
	}
	inline poly& operator*=(const poly &b) { return (*this)=(*this)*b; }
	inline int mul(int *b,int m)
	{
		int t=fast_pow(b[m],mod-2);
		for(int i=n;i>=m;i--)
		{
			int x=(lint)v[i]*t%mod;
			rep(j,0,m) v[i-j]-=(lint)b[m-j]*x%mod,(v[i-j]<0?v[i-j]+=mod:0);
		}
		return n=m-1;
	}
	inline int show()const { debug(n)ln;rep(i,0,n) cerr<<v[i]sp;cerr ln;return 0; }
}ansp;
inline int polymul(char *k,int *b,int n)
{
	poly x;ansp.n=0,ansp.v[0]=1;x.v[0]=0,x.v[1]=1,x.n=1;
	for(int i=(int)strlen(k+1);i;i--,x*=x,x.mul(b,n))
		if(k[i]=='1') ansp*=x,ansp.mul(b,n);return 0;
}
char k[100010];
int main()
{
	scanf("%s",k+1);int n=inn();
	rep(i,1,n) rep(j,1,n) v[i][j]=inn();
	rep(i,0,n)
	{
		b[i]=calc(n,i);
		rep(j,a[i][0]=1,n) a[i][j]=(lint)a[i][j-1]*i%mod;
	}
	gauss(a,b,n),polymul(k,b,n);
	A.init(n),B.init(n),ans.init(n);
	rep(i,1,n) A.v[i][i]=1;
	rep(i,1,n) rep(j,1,n) B.v[i][j]=v[i][j];
	rep(i,0,n-1) ans+=A*ansp.v[i],A*=B;
	return ans.print();
}