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

So Easy!

Time Limit: 1000ms, Special Time Limit:2500ms, Memory Limit:65536KB

Total submit users: 18, Accepted users: 8

Problem 13150 : No special judgement

Problem description

A sequence Sn 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 Sn. 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 < 215, (a-1)2< b < a2, 0 < b, n < 231. The input will finish with the end of file.

Output

For each the case, output an integer Sn.

Sample Input

2 3 1 2013 2 3 2 2013 2 2 1 2013

Sample Output

4 14 4

Problem Source

HNU 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;
}