Number of Battlefields UVA - 11885 矩阵快速幂（高斯消元解线性递推）

Number of Battlefields UVA - 11885

题目大意：给出周长p，问多少种形状的周长为p的，并且该图形的最小包围矩阵的周长也是p，不包括矩形。

假题解：如果包含矩形的话，对应的则是斐波那契数列的偶数项，所以对应减去矩形的个数即可。
讲道理，搜了很多很多博客，都没有讲斐波那契是怎么推的。
可能有人就是想到加上那个包括矩形的数量打个表看看，然后发现了这是斐波那契数列的偶数项。 （我觉得这个可能性大点）

也可能是根据题目给的性质推的。
首先显然奇数项都是0，这里只讨论偶数项。
这里给出一个不是斐波那契的递推，
前面四项：
2 9 29 83
$f_n=5f_{n-1}-8f_{n-2}+5f_{n-3}-f_{n-4}$

它怎么来的呢？
已知前k项求最短递推式可以用Berlekamp-Massey算法
洛谷模板题
模板（看5.3）
所以本题解完结。
.

.

.
但是BM算法好像太难了。
我们通过暴力打表得到了前k项的值$f_k$。如果我们认为它是线性递推，
我们可以设
$f_n=x_1\ast f_{n-1}+x_2\ast f_{n-2}+x_3\ast f_{n-3}+x_4\ast f_{n-4}...$
即有
$$\left(\begin{array}{ccccc}{f}_1& f_2& f_3& ...& f_t\\ f_2& f_3& f_4& ...& f_{t+1}\\ f_3& f_4& f_5& ...& f_{t+2}\\ ...& & & & \\ f_t& f_{t+1}& f_{t+2}& ...& f_{t+t}\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ...\\ x_t\end{array}\right)=\left(\begin{array}{c}{f}_{t+1}\\ f_{t+2}\\ f_{t+3}\\ ...\\ f_{t+t+1}\end{array}\right)$$
所以用高斯消元解如上方程组即可得到$x_1,x_2,x_3,...$的系数，从而求出递推式。
那我们不知道递推式有多少项怎么办。只要我们给出的项足够多。那么上面那个递推方程就会有自由未知元，我们可以令自由未知元为0从而构造出递推方程。
或者直接在高斯消元的过程中，只要找到全0行，那也说明递推式的值已经可以全部确定。（不会证明）
例如，斐波那契递推解的过程：

上图红框处是解的结果。

BM算法求解的时间复杂度是$O(k^2)$，高斯消元显然是$O(k^3)$,$k$是递推式的“长度”。已知递推式求第n项就可以有各种方法啦。

一下是vjudge上ac的代码。

#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#define rep(i,a,b) for(int i=a;i<b;i++)
#define mem(a,b) memset(a,b,sizeof(a))
typedef long long ll;
#include<queue>
using namespace std;
ll mod=1e9+7;
typedef vector<ll> vec;
typedef vector<vec> mat;
ll n,k;
ll qpow(ll x,ll n){//快速幂
	ll c=1;
	for(x%=mod;n;n/=2,x=x*x%mod) if(n&1) c=c*x%mod;
	return c;
}
int N;
void mul(mat&a,mat&b){//矩阵乘法
	mat c(N,vec(N));
	rep(i,0,N) rep(j,0,N) rep(k,0,N)
		(c[i][j]+=a[i][k]*b[k][j])%=mod;
	a=c;
}
void qpow(mat&a,ll n){//矩阵快速幂
	mat gg(N,vec(N));rep(i,0,N) gg[i][i]=1;
	for(;n;n/=2,mul(a,a)) if(n&1) mul(gg,a);
	a=gg;
}
//高斯消元求解线性递推方程
vec gauss(vec xn){//传入的是前多项的系数
	N=xn.size()/2;
	mat dt(N,vec(N));
	rep(i,0,N) rep(j,0,N) dt[i][j]=xn[i+j];//构造矩阵
	int r=0,p=r;
	for(;r<N;r++){
		for(p=r;p<N;p++) if(dt[p][r]) break;//找到该列非0的行
		if(p==N) break;
		swap(dt[p],dt[r]);
		ll ap=qpow(dt[r][r],mod-2);
		rep(j,r,N) dt[r][j]=dt[r][j]*ap%mod;
		rep(i,0,N) if(i!=r){
			ap=dt[i][r];
			rep(j,r,N) dt[i][j]=((dt[i][j]-dt[r][j]*ap)%mod+mod)%mod;
		}
	}
	vec res;
	rep(i,0,r) res.push_back(dt[r-1-i][r]);//把第r列的解拿出来
	return res;
}
//以下是暴力打表求解
int jjj(deque<int> v){
	while(v.size()>1&&v.front()<=v[1]) v.pop_front();//左递增
	while(v.size()>1&&v.back()<=v[v.size()-2]) v.pop_back();//右递减
	return v.size()==1;//单峰
}
ll count_(int x){
	x/=2;
	ll ret=0;
	for(int col=2;col<x-1;col++){//枚举列
		int lin=x-col;
		ll pa=qpow(lin,col);
		for(int i=0;i<pa-1;i++){
			deque<int> d;
			int t=i;
			do d.push_back(t%lin),t/=lin;while(t);
			if(*max_element(d.begin(),d.end())!=lin-1) continue;
			ret+=jjj(d);
		}
	}
	return ret;
}
//以上是暴力打表求解
int main(){//ve是暴力打表结果
	vec ve{2,9,29,83,226,602,1588,4171,10935,28645,75012};
	// for(int i=8;i<30;i+=2) printf("%lld %d\n",count_(i),i),ve.push_back(count_(i));
	vec res=gauss(ve);//求解递推式
	N=res.size();//递推式的长度
	mat init(N,vec(N));//构造递推矩阵
	rep(i,0,N) init[0][i]=res[i]>987654321?(987654321-(mod-res[i])):res[i];
	rep(i,1,N) init[i][i-1]=1;
	mod=987654321;
	ll n;
	while(scanf("%lld",&n),n){
		if(n&1||n<8) {puts("0");continue;}
		n=n/2-3;
		if(n<N) printf("%lld\n",ve[n-1]);
		else {
			mat t=init;
			qpow(t,n-N);
			ll sum=0;
			for(int i=0;i<N;i++) sum=(sum+t[0][i]*ve[N-1-i])%987654321;
			printf("%lld\n",sum);
		}
	}
}

目前这个板子并没有做太多的测试，不能保证正确性。
大概就是懂其原理即可。