HDU3493 The Little Architect 【BM算法推公式+矩阵快速幂】

题目链接：点击打开链接

首先先考虑dp解法。设dp[i][j]为当前使用了i个块，最下一层的块数为j块的时候的方案数。

于是可以得到：

·dp[i][i]=1

·dp[i][j]=dp[i-j][k]*(k-j+1)    (1<=k<=i-j)     (i>j)

于是答案为ans[n]=dp[n][i]   (1<=i<=n)

时间复杂度为O(n^3)，显然超时。

不过可以利用dp的方式打出前几项的表：

1,2,6,19,61,196,629,2017......

有了前几项之后，就可以套用BM算法尝试计算出其递推式：

#include <bits/stdc++.h>

using namespace std;

#define sfi(a) scanf("%d",&a)
#define sfd(a) scanf("%lf",&a)
#define sfl(a) scanf("%lld",&a)
#define sfs(a) scanf("%s",a)

#define rep(i,a,b) for(int i=int(a);i<int(b);++i)
#define dwn(i,b,a) for(int i=int(b-1);i>=int(a);--i)

#define mem(a,p) memset(a,p,sizeof(a))

typedef long long LL;
typedef unsigned UINT;
typedef unsigned long long ULL;

#define MAXN 1005

struct BM
{
	int n;
	
	vector<double> ps[MAXN];
	int pn,fail[MAXN];
	double delta[MAXN];
	
	void Solve(double *x,int n)
	{
		pn=0;
		mem(fail,0);
		mem(delta,0);
		ps[0].clear();
		rep(i,1,n+1)
		{
			double dt=-x[i];
			rep(j,0,ps[pn].size())
				dt+=x[i-j-1]*ps[pn][j];
			delta[i]=dt;
			if(fabs(dt)<=1e-8)continue;
			fail[pn]=i;
			if(!pn)
			{
				ps[++pn].resize(1);
				continue;
			}
			vector<double> &ls=ps[pn-1];
			double k=-dt/delta[fail[pn-1]];
			vector<double> cur;
			cur.resize(i-fail[pn-1]-1);
			cur.push_back(-k);
			rep(j,0,ls.size())cur.push_back(ls[j]*k);
			if(cur.size()<ps[pn].size())cur.resize(ps[pn].size());
			rep(j,0,ps[pn].size())cur[j]+=ps[pn][j];
			ps[++pn]=cur;
		}
	}
	
	void print()
	{
		rep(g,0,ps[pn].size())
			printf("%lf ",ps[pn][g]);
		printf("\n");
	}
}B;

double x[MAXN];

int main()
{
	int n;
	while(sfi(n)==1)
	{
		rep(i,1,n+1)
			sfd(x[i]);
		B.Solve(x,n);
		B.print();
	}
}

发现当输入项数足够多的时候，输出一直为5,-7,4,0

于是便可以得到递推式为ans[n]=5ans[n-1]-7ans[n-2]+4ans[n-3]+0ans[n-4]

矩阵快速幂之即可。

代码：

#include <bits/stdc++.h>

using namespace std;

#define sfi(a) scanf("%d",&a)
#define sfd(a) scanf("%lf",&a)
#define sfl(a) scanf("%lld",&a)
#define sfs(a) scanf("%s",a)

#define rep(i,a,b) for(int i=int(a);i<int(b);++i)
#define dwn(i,b,a) for(int i=int(b-1);i>=int(a);--i)

#define mem(a,p) memset(a,p,sizeof(a))

typedef long long LL;
typedef unsigned UINT;
typedef unsigned long long ULL;

const LL mod=1e9+7;

struct Matrix
{
	LL a[4][4];
	
	void init(LL p)
	{
		mem(a,0);
		rep(i,0,3)a[i][i]=p;
	}
	
	Matrix operator *(const Matrix &p)const
	{
		Matrix ret;
		ret.init(0LL);
		rep(i,0,3)rep(j,0,3)rep(k,0,3)
			ret.a[i][j]=(ret.a[i][j]+a[i][k]*p.a[k][j])%mod;
		return ret;
	}
	
	Matrix operator ^(LL n)
	{
		Matrix k,ret;
		rep(i,0,3)rep(j,0,3)k.a[i][j]=a[i][j];
		ret.init(1LL);
		while(n)
		{
			if(n&1)ret=ret*k;
			k=k*k;
			n>>=1;
		}
		return ret;
	}
};

int main()
{
	LL n;
	while(sfl(n)==1&&n)
	{
		if(n<5)
		{
			if(n==1)printf("1\n");
			if(n==2)printf("2\n");
			if(n==3)printf("6\n");
			if(n==4)printf("19\n");
			continue;
		}
		Matrix M;
		M.init(0LL);
		M.a[0][0]=5LL;
		M.a[0][1]=1e9;
		M.a[0][2]=4LL;
		M.a[1][0]=1LL;
		M.a[2][1]=1LL;
		M=M^(n-4);
		LL ans=M.a[0][0]*19LL+M.a[0][1]*6LL+M.a[0][2]*2LL;
		ans%=mod;
		printf("%lld\n",ans);
	}
}