HDU 4565 So Easy! 递推数列+矩阵快速幂

题目链接点这儿

就是求这个式子的值。。。第一眼。。。太鬼畜了。。。完全是不可做啊。。。

结果。。。练习赛过了会儿之后，，，有人做出来了。。。只好去瞪了。。。发现题目里有个条件(a-1)²< b < a² 也就是说a-1<sqrt(b)<a 这个有什么用呢。。我们可以凑一个跟所求式对称的式子 (a-sqrt(b))^n.。。。这个式子是小于1的，而且这个式子和所求式相加是一个整数（这个凑对称在和二项式相关的问题里还是经常用到的。。）那么Sn就是两个式子的值了。。。

那么，Sn = a^n+b^n。。。这很明显Sn有一个二阶常系数线性齐次的递推公式（一个差分方程），我们将它求出来之后，就可以用矩阵快速幂求第n项了。。

特征根为a+sqrt(b)和a-sqrt(b)。。。两个系数一个是2*a 一个是b-a^2（韦达定理求系数）然后只剩敲一个矩阵乘法了。。

下面放出代码

#include <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <queue>
#include <stack>
#include <bitset>
#include <algorithm>
#include <functional>
#include <numeric>
#include <utility>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <string>
#include <cstring>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <climits>
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)>(b)?(b):(a))
#define rep(i,initial_n,end_n) for(ll (i)=(initial_n);(i)<(end_n);i++)
#define repp(i,initial_n,end_n) for(ll (i)=(initial_n);(i)<=(end_n);(i)++)
#define reep(i,initial_n,end_n) for((i)=(initial_n);(i)<(end_n);i++)
#define reepp(i,initial_n,end_n) for((i)=(initial_n);(i)<=(end_n);(i)++)
#define eps 1.0e-9
#define MAX_N 1010

using namespace std;
typedef pair<int, int> pii;
typedef pair<double, double> pdd;
typedef __int64 ll;
typedef unsigned __int64 ull;


int main(){
    ll a, b, n, m;
    while(scanf("%I64d%I64d%I64d%I64d",&a,&b,&n,&m)!=EOF){
        ll xishu[2][2]={2*a,b-a*a,1,0};
        ll ans[2]={2*a,2};
        while(n) {
            if(n&1) {
                ll a0,a1;
                a0=(xishu[0][0]*ans[0]+xishu[0][1]*ans[1])%m;
                a1=(xishu[1][0]*ans[0]+xishu[1][1]*ans[1])%m;
                ans[0]=a0;
                ans[1]=a1;
            }
            n/=2;
            ll a0,a1,a2,a3;
            a0=(xishu[0][0]*xishu[0][0]+xishu[0][1]*xishu[1][0])%m;
            a1=(xishu[0][0]*xishu[0][1]+xishu[0][1]*xishu[1][1])%m;
            a2=(xishu[1][0]*xishu[0][0]+xishu[1][1]*xishu[1][0])%m;
            a3=(xishu[1][0]*xishu[0][1]+xishu[1][1]*xishu[1][1])%m;
            xishu[0][0]=a0;
            xishu[0][1]=a1;
            xishu[1][0]=a2;
            xishu[1][1]=a3;
        }
        printf("%I64d\n",(ans[1]+m)%m);
    }
    return 0;
}