HDU 1588 Gauss Fibonacci

/*
*  [题意]
*   g(i) = k*i + b
*   f(0)=0, f(1)=1, f(n)=f(n-1)+f(n-2)
*   已知k, b, n, M
*   求( f(g(0))+f(g(1))+...+f(g(n-1)) ) % M
*
*  [解题方法]
*   设斐波那契矩阵A：{1, 1
*                  1, 0}
*   设B = (A^n)
*   则有：B[0][0] = f(n+1), B[0][1] = B[1][0] = f(n), B[1][1] = f(n-1);
*   //上述为斐波那契矩阵的性质
*
*   则只需求：C = (A^b) * ( (A^k)^0+(A^k)^1+(A^k)^2+...+(A^k)^(n-1) ) % M;
*   C[0][1]即为所求
*/
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
using namespace std;
#define LL long long
#define FF(i, n) for(int i = 0; i < n; i++)
#define M 2

int mod;

struct mat {
    int x[M][M];
};

mat matadd(const mat &a, const mat &b)
{
    mat res;
    FF(i, M) FF(j, M)
        res.x[i][j] = (a.x[i][j]+b.x[i][j]) % mod;
    return res;
}

mat matmul(const mat &a, const mat &b)
{
    mat res;
    FF(i, M) FF(j, M) res.x[i][j] = 0;
    FF(i, M) FF(k, M) if(a.x[i][k]) FF(j, M) if(b.x[k][j])
        res.x[i][j] = (res.x[i][j]+(LL)a.x[i][k]*b.x[k][j]%mod) % mod;
    return res;
}

mat qmod(mat a, int b)
{
    mat res;
    FF(i, M) FF(j, M) res.x[i][j] = (i==j);
    for ( ; b; b >>= 1)
    {
        if (b & 1) res = matmul(res, a);
        a = matmul(a, a);
    }
    return res;
}

mat cal(mat a, int n)    //分治求（a^0+a^1+a^2+...+a^n）%mod
{
    if (n == 0) return qmod(a, 0);
    int b = (n+1)/2;
    mat o = cal(a, b-1);
    mat res = matadd(o, matmul(qmod(a, b), o));
    if (n % 2 == 0) res = matadd(res, qmod(a, n));
    return res;
}

int main()
{
    mat A;
    int k, b, n;
    while (cin >> k >> b >> n >> mod)
    {
        //特判
        if (n == 0) {
            puts("0");
            continue;
        }
        //斐波那契矩阵
        A.x[0][0] = A.x[0][1] = A.x[1][0] = 1;
        A.x[1][1] = 0;

        //求(A^b)*((A^k)^0 + (A^k)^1 + ... + (A^k)^(n-1))
        A = matmul(qmod(A, b), cal(qmod(A, k), n-1));

        cout << A.x[0][1] << endl;
    }
    return 0;
}