HDUOJ 4565 So 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<2310< a, m < 2^{15}, (a-1)^2< b < a^2, 0 < b, n < 2^{31}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

非常经典的矩阵快速幂构造题：
首先看题目的条件：
(a−1)2<b<a2(a-1)^2< b < a^2(a−1)2<b<a2
0<a−b<10<a-\sqrt{b}<10<a−b<1 => 可以推出 a,ba,ba,b 必有一个是小数
(a−b)n<1(a-\sqrt{b})^n<1(a−b)n<1
┌(a−b)n(a-\sqrt{b})^n(a−b)n┐=1=1=1 => ┌(a+b)n(a+\sqrt{b})^n(a+b)n┐=(a+b)n+(a−b)n(a+\sqrt{b})^n+(a-\sqrt{b})^n(a+b)n+(a−b)n
所以我们考虑构造:
s=(a+b)n+(a−b)ns=(a+\sqrt{b})^n+(a-\sqrt{b})^ns=(a+b)n+(a−b)n
下面找递推矩阵：
类比斐波那契数列，我们要找 Sn=pSn−1+qSn−2S_n=pS_{n-1}+qS_{n-2}Sn=pSn−1+qSn−2
移项得推导方程：
S2−pS−q=0S^2-pS-q=0S2−pS−q=0
已知两解为：
S1=a+b,S2=a−bS_1=a+\sqrt{b},S_2=a-\sqrt{b}S1=a+b,S2=a−b
解得：
p=2a,q=a−b2p=2a,q=a-b^2p=2a,q=a−b2
所以得到递推矩阵：
[SnSn−1]=[2aa−b210][Sn−1Sn−2]\left[ \begin{matrix} S_n\\ S_{n-1} \end{matrix} \right]= \left[ \begin{matrix} 2a&a-b^2\\ 1&0\\ \end{matrix} \right] \left[ \begin{matrix} S_{n-1}\\ S_{n-2} \end{matrix} \right] [SnSn−1]=[2a1a−b20][Sn−1Sn−2]
AC代码如下：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll mod,n,A,B;
struct mat
{ll m[2][2];
}ans,a;void init(mat &a){memset(a.m,0,sizeof(a.m));for(int i=0;i<2;i++) a.m[i][i]=1;
}mat mul(mat a,mat b)
{mat c;memset(c.m,0,sizeof(c.m));for(int i=0;i<2;i++){for(int k=0;k<2;k++){if(a.m[i][k]==0) continue;for(int j=0;j<2;j++){c.m[i][j]=(c.m[i][j]+a.m[i][k]*b.m[k][j])%mod;if(c.m[i][j]<0) c.m[i][j]+=mod;}}}return c;
}mat mat_pow(mat a,ll k)
{mat ans;init(ans);while(k>0){if(k&1) ans=mul(a,ans);a=mul(a,a);k>>=1;}return ans;
}int main(){while(cin>>A>>B>>n>>mod){a.m[0][0]=2*A%mod;a.m[0][1]=((B-A*A%mod)+mod)%mod;a.m[1][0]=1;a.m[1][1]=0;ans.m[0][0]=2*A%mod;ans.m[1][0]=2%mod;ans=mul(mat_pow(a,n-1),ans);cout<<ans.m[0][0]<<endl;}return 0;
}