hdu 1588(Fibonacci矩阵求和)

题目的大意就是求等差数列对应的Fibonacci数值的和,容易知道Fibonacci对应的矩阵为[1,1,1,0],因为题目中f[0]=0,f[1]=1,所以推出最后结果f[n]=(A^n-1).a,所以 f(g(i))= f(k*i+b)= (A^(k*i+b-1)).a,i从 0取到 n-1,取出公因式 A^(b-1)(因为矩阵满足分配率),然后所求结果可化为 A^(b-1) * (A^0 + A^k + A^2k +....+ A^(n-1)k),化到这里后难点就是求和了,一开始我尝试暴力求和(每个A^k可以用快速幂求出,logn级别),即O(n)的做法,结果TLE了,预料之中,这时我竟傻乎乎地套用等比数列求和公式,即(A^nk -A^0) /(A^k -E),按比例放大再相减是没错,问题是不是简单的相除……总之思路应该是错的了,后来看别人的博客后才知道原来可以用二分来求和的,即 A^0 + A^k + A^2k +....+ A^(n-1)k = (A^0 + A^k + A^2k +...+ A^(n/2-1)k) *(E + A^(n/2)k),递归来求和,处理好 n的奇偶性即可,但我还是调试了好久,就因为在一些细节问题上出错却检查不出来,下面附上代码


#include<cstdio>
#include<cstring>
#include<cstdlib>
using namespace std;
typedef long long LL;
LL mod= 1000000007;		//初不初始化都没问题,只是为了防止忘记读入时产生的异常退出 

struct matrix{
	LL a,b,c,d;
	matrix(LL a=0, LL b=0, LL c=0, LL d=0): a(a),b(b),c(c),d(d) {}
	matrix operator +(const matrix &m2){
		return matrix((a+m2.a)%mod,(b+m2.b)%mod,(c+m2.c)%mod,(d+m2.d)%mod);
	}
	matrix operator *(const matrix &m2){
		return matrix((a*m2.a%mod+b*m2.c%mod)%mod, (a*m2.b%mod+b*m2.d%mod)%mod, (c*m2.a%mod+d*m2.c%mod)%mod, (c*m2.b%mod+d*m2.d%mod)%mod);
	}
	//一开始等比数列求和的思路要用到除法,后来才发现是错的,不过也不删了,就放在这吧 
	matrix operator /(const matrix &m2){
		//二维求逆很好求的 
		matrix inv= matrix(m2.d,-m2.b,-m2.c,m2.a);
		LL tmp= m2.a*m2.d-m2.b*m2.c;
		return matrix((a*inv.a/tmp%mod+b*inv.c/tmp%mod)%mod, (a*inv.b/tmp%mod+b*inv.d/tmp%mod)%mod, (c*inv.a/tmp%mod+d*inv.c/tmp%mod)%mod, (c*inv.b/tmp%mod+d*inv.d/tmp%mod)%mod);
	}
}; 
// A为 Fibonacci矩阵,E为单位矩阵,设为全局变量更方便一些 
matrix A(1,1,1,0),E(1,0,0,1);

//简单的快速幂 
matrix quick_mod(matrix m, LL b){
	if(b==-1)	return matrix(0,1,1,-1);
	//若指数为-1,返回矩阵 A^-1,相当于A^1的逆(算了好久T.T)
	matrix res(E);		//res一开始为单位矩阵 
	while(b){
		if(b&1)	  res= res*m;
		m= m*m;
		b>>=1;
	}
	return res;
}

//二分法计算 A^0 + A^k + A^2k +....+ A^(n-1)k 的和
//即 sum =(A^0 + A^k + A^2k +...+ A^(n/2-1)k) *(E + A^(n/2)k),分治的思想 
matrix quick_sum(LL k, LL n){
	if(n==1)	return E;	//若 n为1,即计算 A^0,此时返回的是 E!而不是 A!一开始没想到错在这里,卡了好久 T.T 
	if(n%2==0)	return quick_sum(k,n/2)*(E+quick_mod(A,(n/2)*k));
	//else	return quick_sum(k,(n-1)/2)*(E+quick_mod(A,((n-1)/2)*k))*quick_mod(A,(n-1)*k);
	else	return quick_sum(k,n-1) + quick_mod(A,(n-1)*k);
	//这样的写法比起上一行的虽然会多一次调用函数的开销,但可读性增强,代码逻辑更清晰 
}

int main()
{
	//freopen("1588in.txt","r",stdin);
	LL k,b,n;
	while(~scanf("%I64d%I64d%I64d%I64d",&k,&b,&n,&mod)){
		//最后的结果就是 A^b-1 *(A0 + A^k + A^2k +...+ A(n-1)k) 的 a
		matrix ans= quick_mod(A,b-1) * quick_sum(k,n);
		printf("%I64d\n",ans.a);
	}
	return 0;
}

  • 1
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值