HDU 1588 Gauss Fibonacci 矩陣

[cpp] 
/**
  *stratege: 構造矩陣 A =  |1 1|
  *                        |1 0|
  *          矩陣快速冪，二分求和
  *status: johnsondu B Accepted 252 KB 15 ms C++ 3113 B 
  *
*/ 
 
 
#include <iostream> 
#include <cstdio> 
#include <cstring> 
#include <cmath> 
#include <algorithm> 
 
using namespace std; 
 
typedef __int64 LLD ; 
 
struct Matrix  
{ 
    LLD mat[2][2] ; 
}; 
Matrix tmp, unit, res ;  
LLD K, B, n, M ; 
Matrix a1, a2, a3 ; 
 
Matrix Multi (Matrix a, Matrix b)  //兩矩陣相乘 
{ 
    Matrix c ; 
    for (int i = 0; i < 2;i ++) 
        for (int j = 0; j < 2; j ++) 
        { 
            c.mat[i][j] = 0 ; 
            for (int k = 0; k < 2; k ++) 
                c.mat[i][j] += (LLD)a.mat[i][k]*b.mat[k][j] ,c.mat[i][j] %= M ; //矩陣相乘時，需要防止溢出 
        } 
    return c ; 
} 
 
Matrix power (LLD p)  //求出a1 = A^b 
{ 
    Matrix tp = unit ; 
    Matrix tr = tmp ; 
    while (p) 
    { 
        if (p&1) 
        { 
            tr = Multi (tr, tp) ; 
            p -- ; 
            continue ; 
        } 
        tp = Multi (tp, tp) ; 
        p >>= 1 ; 
    } 
    return tr ; 
} 
 
 
Matrix power1 (LLD p) //求出a2 = A^k 
{ 
    Matrix tp = a2 ; 
    Matrix tr = tmp ; 
    while (p) 
    { 
        if (p&1) 
        { 
            tr = Multi (tr, tp) ; 
            p -- ; 
            continue ; 
        } 
        tp = Multi (tp, tp) ; 
        p >>= 1 ; 
    } 
    return tr ; 
} 
 
Matrix add (Matrix a, Matrix b)  //矩陣相加 
{ 
    Matrix c ; 
     
    for (int i = 0; i < 2; i ++) 
        for (int j = 0; j < 2; j ++) 
            c.mat[i][j] = (a.mat[i][j] + b.mat[i][j]) % M ; 
    return c ; 
} 
 
void init () 
{ 
    unit.mat[0][0] = 1 ;  //初始化矩陣 
    unit.mat[0][1] = 1 ; 
    unit.mat[1][0] = 1 ; 
    unit.mat[1][1] = 0 ; 
    tmp.mat[0][0] = 1 ;   //單位矩陣 
    tmp.mat[0][1] = 0 ; 
    tmp.mat[1][0] = 0 ; 
    tmp.mat[1][1] = 1 ; 
    a1 = power (B) ;      //由題意有ans = A^b + A^(k+b) + A^(2k+b) + ... + A^((n-1)k + b) 
    a2 = power (K) ;      //由此得出 a1 = A^b, a2 = A^k, 所以有ans = a1 * (E + a2^1 + a2^2 + ... + a2^(n-1)) ; 
} 
 
Matrix MatrixSum (LLD p)  //求出E + a2^1 + a2^2 + ... + a2^(n-1) 
{ 
    if (p == 1)           //注意，此時p==1時，基為a2 = A^k 
        return a2 ;  
    Matrix tm, tr ; 
    tm = MatrixSum (p/2) ; 
    if (p & 1) 
    { 
        tr = power1 (p/2+1) ; 
        tm = add (tm, Multi (tm, tr)) ; 
        tm = add (tm, tr) ; 
    } 
    else { 
        tr = power1 (p/2) ; 
        tm = add (tm, Multi (tm, tr)) ; 
    } 
    return tm ; 
     
} 
 
void output ()  //輸出 
{ 
    if (n == 1) 
    { 
        printf ("%I64d\n", a1.mat[0][1]) ; 
        return ; 
    } 
         
    if (n == 0) 
    { 
        printf ("0\n") ; 
        return ; 
    } 
     
    Matrix ans = MatrixSum (n-1) ;  //因為是求出0~n-1, 此處先求出1~n-1 
    ans = add (ans, tmp) ;          //再求出與單位矩陣的和 
    ans = Multi (ans, a1) ;         //ans = a1 * (E + a2^1 + a2^2 + ... + a2^(n-1)) 
    printf ("%I64d\n", ans.mat[0][1] % M) ; //輸出 
} 
 
int main() 
{ 
    while (scanf ("%I64d%I64d%I64d%I64d", &K, &B, &n, &M) != EOF) 
    { 
        init () ; 
        output () ; 
    } 
    return 0; 
}