hnu13150（hdu 4565）SO EASY!（矩阵快速幂）

来源：互联网发布：vb.net 字符串包含编辑：程序博客网时间：2024/05/29 16:05

So Easy!Time Limit: 1000ms, Special Time Limit:2500ms, Memory Limit:65536KBTotal submit users: 18, Accepted users: 8Problem 13150 : No special judgementProblem description

A sequence S_n is defined as:

Where a, b, n, m are positive integers.┌x┐is the ceil of x. For example, ┌3.14┐=4. You are to calculate S_n. You, a top coder, say: So easy!

Input

There are several test cases, each test case in one line contains four positive integers: a, b, n, m. Where 0< a, m < 2¹⁵, (a-1)²< b < a², 0 < b, n < 2³¹.
The input will finish with the end of file.

Output

For each the case, output an integer S_n.

Sample Input

2 3 1 20132 3 2 20132 2 1 2013

Sample Output

Problem SourceHNU Contest

原来博士出的SO EASY就是这个题。。。之前听过了但是没做过

之前做过一个类似的所以一看到就知道是用矩阵快速幂做

但是这个多了一个取上界第一直觉就是按之前的想法求值之后加一

按照这个想法写了之后跑了样例一样。。。就交了发现WA了

就意识到应该是那个上界的处理上出了问题

因为 ceil((a+sqrt(b))^n)=(a+sqrt(b))^n+(a-sqrt(b)^n)

令 a+sqrt(b)=x a-sqrt(b)=y

则(a+sqrt(b))^n+(a-sqrt(b)^n)=x^n+y^n=(x+y)(x^(n-1)+y^(n-1))-x*y(x^(n-2)+y^(n-2))

设 g(n)=x^n+y^n x+y=2*a=p x*y=a*a-b=q

|g(n)g(n-1)|=|g(n-1) g(n-2)|* |p 1|

|-q 0|

然后就可以由g(1) g(0) 2阶矩阵的n-1次幂求值

g(1)=2*a=p g(0)=2

注意数据范围所有的都要是__int64

代码如下：

#include <cstdio>#include <iostream>#include <cstring>#include <cmath>#include <algorithm>#include <string.h>#include <string>#define MEM(a,x) memset(a,x,sizeof a)#define eps 1e-8//#define MOD 10009#define MAXN 10010#define INF 0x7fffffff#define ll __int64#define bug cout<<"here"<<endl;#define fread freopen("ceshi.txt","r",stdin)#define fwrite freopen("out.txt","w",stdout)using namespace std;ll mod,c,d;struct matrix{    ll a[2][2];};matrix ori,res;void init(){    res.a[0][1]=res.a[1][0]=0;    res.a[0][0]=res.a[1][1]=1;}matrix multiply(matrix x,matrix y){    matrix temp;    memset(temp.a,0,sizeof(temp.a));    for(int i=0;i<2;i++)        for(int j=0;j<2;j++)            for(int k=0;k<2;k++)    {        temp.a[i][j]+=x.a[i][k]*y.a[k][j];        temp.a[i][j]%=mod;    }    return temp;}void calc(ll n){    while(n)    {        if(n&1)   res=multiply(ori,res);        n>>=1;        ori=multiply(ori,ori);    }}int main(){//    fread;    ll n;    while(scanf("%I64d%I64d%I64d%I64d",&c,&d,&n,&mod)!=EOF)    {        init();        ll x=2*c;        ll y=c*c-d;        ori.a[0][0]=x;        ori.a[1][1]=0;        ori.a[0][1]=1;        ori.a[1][0]=-y;//        ori.a[0][0]=ori.a[1][1]=c;//        ori.a[0][1]=d;//        ori.a[1][0]=1;        calc(n-1);        ll ans;        ans=((x*res.a[0][0]+2*res.a[1][0])%mod+mod)%mod;//        double ans=0.0;//        ll an=res.a[0][0]*c%mod+res.a[0][1]%mod;//        an%=mod;//        ll bn=res.a[1][0]*c%mod+res.a[1][1]%mod;//        bn%=mod;//        ans+=(double)1.0*res.a[0][0]*c%mod;//        ans+=(double)1.0*res.a[0][1]%mod;//        ans+=(double)(res.a[1][0]*c*1.0+1.0*res.a[1][1])*sqrt((double)5)%mod;//        ans=ceil((double)((double)an+(double)bn*1.0*sqrt((double)d)));//        printf("%llf\n",ans);//        ans=(double)an+(double)bn*1.0*sqrt((double)d);//        ll pp=(int)ans;//        pp++;//        pp%=mod;        printf("%I64d\n",ans);    }    return 0;}

0 0