題目鏈接：http://acm.hdu.edu.cn/showproblem.php?pid=4565

題解：(a+√b)^n=xn+yn*√b,(a-√b)^n=xn-yn*√b,

(a+√b)^n=2*xn-(a-√b)^n,(0<=a-√b<=1),所以整數部分就是2*xn

然後再利用兩個公式

(a+√b)^(n+1)=(a+√b)*(xn+yn*√b)

(a-√b)^(n+1)=(a-√b)*(xn-yn*√b)

聯立得到

x(n+1)=a*xn+b*yn

y(n+1)=xn+a*yn;

然後就是矩陣快速冪

#include <iostream>
#include <cstring>
#include <cstdio>
using namespace std;
typedef long long ll;
struct Matrix {
    ll dp[3][3];
};
ll a , b , n , m;
Matrix Mul(Matrix a , Matrix b) {
    Matrix c;
    memset(c.dp , 0 , sizeof(c.dp));
    for(int i = 0 ; i < 2 ; i++) {
        for(int j = 0 ; j < 2 ; j++) {
            for(int k = 0 ; k < 2 ; k++) {
                c.dp[i][j] += ((a.dp[i][k] * b.dp[k][j]) % m + m) % m;
            }
        }
    }
    return c;
}
Matrix quick_pow(Matrix a , ll k) {
    Matrix res;
    memset(res.dp , 0 , sizeof(res.dp));
    for(int i = 0 ; i < 2 ; i++) {
        res.dp[i][i] = 1;
    }
    while(k) {
        if(k & 1) res = Mul(res , a);
        k >>= 1;
        a = Mul(a , a);
    }
    return res;
}
int main() {
    while(~scanf("%lld%lld%lld%lld" , &a , &b , &n , &m)) {
        if(n == 0) {
            printf("%lld\n" , (ll)1 % m);
        }
        else {
            Matrix cnt , ans , sta;
            sta.dp[0][0] = 1 , sta.dp[1][0] = 0;
            cnt.dp[0][0] = a , cnt.dp[0][1] = b;
            cnt.dp[1][0] = 1 , cnt.dp[1][1] = a;
            ans = quick_pow(cnt , n);
            ans = Mul(ans , sta);
            printf("%lld\n" , 2 * ans.dp[0][0] % m);
        }
    }
    return 0;
}

hdu4565 So Easy!（矩陣快速冪）