51nod 1358 浮波那契 (类斐波那契数列+矩阵快速幂+构造矩阵)

传送门：51nod 1358

题目大意：

一开始以为 n 只可能是整数……所以看到 FB( n-3.4 ) 的时候一 lemon 逼。这个就是给出一个类斐波那契数列，让你求某一项的值，不同的是它的递推公式里面有小数。

参考资料：Matrix神犇讲解：十个利用矩阵乘法解决的经典题目。

前置技能：

看到这道题不知道大家有没有想到一道用矩阵快速幂求斐波那契数列第 n 项值的题。建议先做下这道题，附上题解：矩阵快速幂求斐波那契数列第 n 项值题解。

方法简单的说就是构造一个方阵，然后让这个方阵每自乘一次后就会让数列的前 2 项的值迭代一次。最后乘以初始化矩阵（即数列前两项），然后结果矩阵中所有数之和就是答案。

所以，我们可不可以按照这个思路来解决这个问题呢？当然是可以的，不然我说它作甚……既然如此，关键就在于如何构造这个方阵和初始矩阵。下面主要讲解这个问题，我的表述和网上有些人的表述有些不同，但是作用是一样的，并且他们在实现的时候也是按照的以下思路。请大家耐心读下去。

首先我们规定，行数从上到下递增，列数从左到右递增。

若类斐波那契数列形式为：

f(N) = a1*f(N-x1) + a2*f(N-x2) + … + an*f(N-xn) , ( x1<x2<…<xn )

f(N) = b1 ,b2 ,… ,bi ,… ,bn , ( 1<= i <=xn )

则我们需要构造一个 n * n 的方阵，其中左下角 (n-1)*(n-1) 的方阵的主对角线为 1 。第 1 行对应的  xi ( 1<= i <=n ）处值为 ai ,其他地方为 0。求第 i (i>xn) 项的值的时候，只需要将构造的方阵自乘 i-xn 次，然后乘以初始矩阵后结果矩阵的第 (1,1) 位就是第 i 项值了。而初始矩阵是根据 1<= i <=xn 时候的取值来的。

例如有 f(N) = 4f(N-1) - 3f(N-2) + 2f(N-4)
            f(N) = 1 ,2 ,3 ,4  ( N=1 ,2 ,3 ,4 )

则其对应的方阵为：             4  -3   0   2
                                             1   0   0   0
                                             0   1   0   0
                                             0   0   1   0

初始矩阵为：        4
                             3
                             2
                             1

求第 i (i>4) 项的时候，将以上方阵自乘 i-4 次，然后乘以初始矩阵，结果矩阵必定是只有 1 列，结果数列的第 (1,1) 位就是第 i 项的值。

思路：

好了，该掌握的知识说的差不多了，下面回到本题。首先要解决的是小数的问题。我们可以构造一个新的函数。让其为：   G(n) = 1 , n<=40

          G(n) = G(n-10) + G(n-34) , n>40

很明显 G(n) = FB( n*10 ) .夹克老爷也说了，其实还可以更优，用 5 和 17 会更好，其实也就是将 40，10，34 同时除以他们的最大公约数 2，具体的就不说了，留给大家思考。下面还是以 10 和 34 来讲解。

下面就是构造方阵和初始化矩阵了，按照上面说的，应该是为 34*34 的方阵，左下角 33*33 的方阵的主对角线为 1 ，第 1 行第 10 列和第 34 列为 1，其他地方为 0. 初始化矩阵很明显为 34*1 的矩阵，全为 1.

再然后如果是求前 40 项（针对于 G(n)来说）就直接输出 1，其他的第 i 项，先将方阵自乘 i-40 次，理论上需要乘以初始矩阵，不过鉴于初始矩阵都是 1 ，所以可以将自乘后矩阵的第 1 行数相加就是答案了。

代码：

#include<stdio.h>
#include<string.h>
typedef long long LL;

LL m=1e9+7;

int n=34; //方阵的行数 

struct matrix
{ //矩阵结构体 
	LL map[40][40];
}a,per;

void init()
{ //初始化矩阵 
	int i,j;
	memset(a.map,0,sizeof(a.map));
	a.map[1][10]=1;  //第 1 行第10和34列为1 
	a.map[1][34]=1;
	for(i=2;i<=n;i++) //左下角小矩阵主对角线为 1 
		a.map[i][i-1]=1;
	for(i=1;i<=n;i++)
		for(j=1;j<=n;j++)
			per.map[i][j]=(i==j); //主对角线为1，其他为0，为单位矩阵
}

matrix multi(matrix a,matrix b)
{ //实现a、b两个矩阵的相乘 
	matrix c;
	int i,j,k;
	for(i=1;i<=n;i++)
		for(j=1;j<=n;j++)
		{
			c.map[i][j]=0;
			for(k=1;k<=n;k++) c.map[i][j]+=(a.map[i][k]*b.map[k][j])%m;
			c.map[i][j]%=m;
		}
	return c;
}

//原理和整数快速幂一样，只是将整数乘法改为了矩阵乘法 
matrix power(LL k)
{ //计算a矩阵的k次幂 
	matrix ans,base;
	ans=per; //初始化为单位矩阵 
	base=a;
	while(k)
	{
		if(k&1) ans=multi(ans,base);
		base=multi(base,base);
		k>>=1;
	}
	return ans;
}

int main()
{
	LL i,t,ans;
	matrix aa;
	init();
	while(~scanf("%lld",&t))
	{		
		t=t*10; //这里的 t 对应的是函数 G(n) 
		if(t<=40)
		{ //如果 <=40 项直接输出 
			printf("1\n");
			continue; 
		}
		t-=40; //不要忘记先减去40 
		aa=power(t); //矩阵快速幂 
		ans=0;
		for(i=1;i<=n;i++)
		{ //将第1行数相加，不要忘记同时取模 
			ans+=aa.map[1][i]%m;
			ans%=m;
		}
		printf("%lld\n",ans);
	}
	return 0;
}