注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

fudq's AC Road

何以解忧,唯有AC!

 
 
 

日志

 
 

hdu 1588 Gauss Fibonacci  

2013-12-23 14:53:02|  分类: ACM-hdu |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
http://acm.hdu.edu.cn/showproblem.php?pid=1588
题意:
计算sigma(f(g(i))),0<=i<n,
g(i)=k*i+b, f(i)为第i项Fibonacci数。
先已知k,b,n,m,求结果对m取余的结果。
题解:
矩阵构造+矩阵快速幂。
展开下得到 f(b)+f(k+b)+f(2*k+b)+f(3*k+b)+……+f((n-1)*k+b)
第i项Fibonacci可以由矩阵快速幂求解得到,
构造矩阵A=| 0 1 |,令f(i)=A^i,f(i)的结果为矩阵A[0, 0]项,化简得A^b+A^b*A^k+A^b*A^2k+A^b*A^3k+……+A^b*A^(n-1)k
                 | 1 1 |
提取A^b,得到A^b(E+A^k+A^2k+A^3k+……+A^(n-1)k),E表示单位矩阵
令T=A^k,得到A^b(E+T+T^2+T^3+……+T^(n-1)),g(n)=E+T+T^2+T^3+……+T^(n-1)
可以发现g(n)可以由g(n-1)得到,g(n)=g(n-1)*T+E,所以可以构造矩阵B,
矩阵B=| E E |,令g(n)=B^n,g(n)的结果为矩阵B[0, 1]项。
     | 0  T |
最后结果Ans=A^b * B^n,这两部分都可以由矩阵快速幂求得。
因为Fibonacci数第一项为0,第二项为1,(0, 1)*Ans的第一项为最后的结果,也就是Ans矩阵的[1 ,0]项。

LL k,b,n,m;

struct Matrix{
LL v[5][5];
};

Matrix Mat_mul(Matrix m1,Matrix m2,LL pri,LL len) //矩阵相乘
{
Matrix c;
MEM(c.v);
for(int i=0;i<len;i++)
for(int j=0;j<len;j++)
for(int k=0;k<len;k++)
c.v[i][j]=(c.v[i][j]+m1.v[i][k]*m2.v[k][j])%pri;
return c;
}

Matrix Mpow(Matrix A,LL n,LL pri,LL len) //矩阵快速幂
{
Matrix c,x=A;
MEM(c.v);
for(int i=0;i<len;i++)
c.v[i][i]=1;
while(n >= 1)
{
if(n & 1)
c=Mat_mul(c,x,pri,len);
n>>=1;
x=Mat_mul(x,x,pri,len);
}
return c;
}

//定义01矩阵tmp
void init(Matrix &tmp)
{
MEM(tmp.v);
tmp.v[0][0]=0;tmp.v[0][1]=1;tmp.v[1][0]=1;tmp.v[1][1]=1;
}

void solve()
{
Matrix ab,bn,tmp,tmp2,ans1,ans2,res;
init(ab);
ans1=Mpow(ab,b,m,2); //A^b
tmp=Mpow(ab,k,m,2); //A^k

//make matrix bn
bn.v[0][0]=1;bn.v[0][1]=0;bn.v[0][2]=1;bn.v[0][3]=0;
bn.v[1][0]=0;bn.v[1][1]=1;bn.v[1][2]=0;bn.v[1][3]=1;
bn.v[2][0]=0;bn.v[2][1]=0;bn.v[2][2]=tmp.v[0][0];bn.v[2][3]=tmp.v[0][1];
bn.v[3][0]=0;bn.v[3][1]=0;bn.v[3][2]=tmp.v[1][0];bn.v[3][3]=tmp.v[1][1];

ans2=Mpow(bn,n,m,4); //B^n
tmp2.v[0][0]=ans2.v[0][2];tmp2.v[0][1]=ans2.v[0][3];tmp2.v[1][0]=ans2.v[1][2];tmp2.v[1][1]=ans2.v[1][3];
res=Mat_mul(ans1,tmp2,m,2);//A^b * B^n
pf64(res.v[1][0]);
}

int main()
{
#ifndef ONLINE_JUDGE
freopen("testin.txt", "r", stdin);
// freopen("textout.txt", "w", stdout);
#endif
while(scanf("%I64d%I64d%I64d%I64d",&k,&b,&n,&m)!=EOF)
solve();
return 0;
}


  评论这张
 
阅读(52)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017