# hdu4565So Easy!+矩阵快速幂+广义斐波拉契

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

Source
2013 ACM-ICPC长沙赛区全国邀请赛——题目重现

(a1)2<b<a2$(a-1)^2< b < a^2$所以依然可以(a+b)n$(a+\sqrt b)^n$(ab)n$(a-\sqrt b)^n$相加。

#include<cstdio>
#include<cmath>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<queue>
#include<vector>
#include<map>
#include<stack>
#include<set>
using namespace std;
#define pi acos(-1.0)
#define EPS 1e-6    //log(x)
#define e exp(1.0); //2.718281828
//#define mod 1000000007
#define INF 0x7fffffff
#define inf 0x3f3f3f3f
typedef long long LL;

#define debug(x) cout<<x<<endl;
#define debug2(x) cout<<x<<" ";

//#define MOD 10000007
LL MOD;
struct Mat{
int n,m;
LL mat[9][9];
};
Mat operator *(Mat a,Mat b){
Mat c;
memset(c.mat,0,sizeof(c.mat));
c.n = a.n,c.m = b.m;

for(int i=1;i<=a.n;i++){
for(int j=1;j<=b.m;j++){
for(int k=1;k<=a.m;k++){
c.mat[i][j] += (a.mat[i][k]*b.mat[k][j])%MOD;
c.mat[i][j] %= MOD;
}
}
}
return c;
}
Mat operator +(Mat a,Mat b){
Mat c;
memset(c.mat,0,sizeof(c.mat));
c.n = a.n,c.m = a.m;

for(int i=1;i<=a.n;i++){
for(int j=1;j<=a.m;j++){
c.mat[i][j] = (a.mat[i][j]+b.mat[i][j])%MOD;
}
}
return c;
}
Mat operator ^(Mat a,LL k){
Mat c;
memset(c.mat,0,sizeof(c.mat));
c.n = a.n,c.m = a.n;
for(int i=1;i<=a.n;i++)c.mat[i][i] = 1;

while(k){
if(k&1){
c = c*a;
}
a = a*a;
k>>=1;
}
return c;
}
void out(Mat a){
for(int i=1;i<=a.n;i++){
for(int j=1;j<=a.m;j++){
printf(j==a.m? "%I64d\n":"%I64d ",a.mat[i][j]);
}
}
}
LL quickPow(LL x, LL n, LL mm)
{
LL a = 1;
while (n)
{
a *= n&1 ? x : 1;
a %= mm;
n >>= 1 ;
x *= x;
x %= mm;
}
return a;
}
int main()
{
LL a,b,n;
while(scanf("%I64d %I64d %I64d %I64d",&a,&b,&n,&MOD)!=EOF){
if(n==1){
printf("%I64d\n",((LL)(a+ceil(sqrt(b*1.0)))%MOD));
continue;
}
Mat pp;
pp.n=pp.m=2;
pp.mat[1][1]=a%MOD;
pp.mat[1][2]=b%MOD;
pp.mat[2][1]=1;
pp.mat[2][2]=a%MOD;

Mat A0;
A0.n=2,A0.m=1;
A0.mat[1][1]=a%MOD;
A0.mat[2][1]=1;

Mat ans=pp^(n-1);
ans=ans*A0;
printf("%I64d\n",(2*ans.mat[1][1])%MOD);
}
return 0;
}

/*
_ooOoo_
o8888888o
88" . "88
(| -_- |)
O\  =  /O
____/---'\____
.'  \\|     |//  .
/  \\|||  :  |||//  \
/  _||||| -:- |||||-  \
|   | \\\  -  /// |   |
| \_|  ''\---/''  |   |
\  .-\__  -  ___/-. /
___. .'  /--.--\  . . __
."" '<  .___\_<|>_/___.'  >'"".
| | :  - \.;\ _ /;./ -  : | |
\  \ -.   \_ __\ /__ _/   .- /  /
======-.____-.___\_____/___.-____.-'======
=---='
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
I have a dream!A AC deram!!
orz orz orz orz orz orz orz orz orz orz orz
orz orz orz orz orz orz orz orz orz orz orz
orz orz orz orz orz orz orz orz orz orz orz

*/

`

#### hdu5451（矩阵快速幂+广义斐波拉契）

2017-12-10 17:06:30

#### 算法学习笔记（五） 递归之 快速幂、斐波那契矩阵加速

2014-08-17 20:19:28

#### 矩阵快速幂求斐波那契数列（初学整理）

2016-10-31 13:27:45

#### hdu2256Problem of Precision+矩阵快速幂+广义斐波拉契

2016-08-31 10:53:36

#### 快速幂与快速矩阵幂(以大数下的斐波那契数列为例)

2017-08-30 22:03:07

#### ,快速乘,快速幂,矩阵快速幂(求斐波那契数列)

2016-12-10 00:19:24

#### poj 3070 Fibonacci 【矩阵快速幂 求第N个斐波那契数%1000】

2015-08-30 20:01:11

#### 斐波那契数列 打表+矩阵快速幂

2016-10-01 20:26:22

#### 算法——矩阵快速幂 求第N个斐波那契数

2017-10-10 20:47:09

#### HDU-4565 So Easy! 矩阵快速幂 & 共轭构造

2017-08-15 00:59:25