HDU 1588 Gauss Fibonacci

HDU_1588

    不妨设F为像下图一样的含有fibonacci项的矩阵，那么最后结果就是矩阵S(n-1)=F^(b)+F^(k+b)+…+F^(k*(n-1)+b)左上角的元素。

    对于矩阵S和F，可以得到递推公式S(n)=S(n-1)+F^(k*n+b)，F^(k*(n+1)+b)=A^(k)*F^(k*n+b)，矩阵A为下图所示的矩阵。于是我们就可以将这个两个递推式构造成一个递推关系的矩阵，然后再用二分矩阵的方法快速求得S(n-1)即可。

    

#include<stdio.h>
#include<string.h>
#define MAXD 8
int K, B, N, M;
struct Matrix
{
    int a[MAXD][MAXD];
    Matrix()
    {
        memset(a, 0, sizeof(a));
    }
};
Matrix multiply(Matrix &x, Matrix &y)
{
    int i, j, k;
    Matrix z;
    for(k = 0; k < 8; k ++)
        for(i = 0; i < 8; i ++)
            if(x.a[i][k])
            {
                for(j = 0; j < 8; j ++)
                    if(y.a[k][j])
                        z.a[i][j] = (z.a[i][j] + (long long)x.a[i][k] * y.a[k][j]) % M;
            }
    return z;
}
Matrix powmod(Matrix &unit, Matrix &mat, int n)
{
    while(n)
    {
        if(n & 1)
            unit = multiply(mat, unit);
        n >>= 1;
        mat = multiply(mat, mat);
    }
    return unit;
}
Matrix getF(int n)
{
    Matrix unit, mat;
    unit.a[0][0] = 1;
    mat.a[0][0] = mat.a[0][1] = mat.a[1][0] = 1;
    powmod(unit, mat, n - 1);
    return unit;
}
Matrix getA(int n)
{
    Matrix unit, mat;
    unit.a[0][0] = unit.a[1][1] = 1;
    mat.a[0][0] = mat.a[0][1] = mat.a[1][0] = 1;
    unit = powmod(unit, mat, n);
    return unit;
}
void solve()
{
    int i, j;
    Matrix unit, mat, Ak, Ab, Akb;
    if(B)
    {
        Ab = getF(B);
        for(i = 0; i < 4; i ++)
            for(j = 0; j < 4; j ++)
                unit.a[i][j] = Ab.a[i][j];
    }
    if(K + B)
    {
        Akb = getF(B + K);
        for(i = 0; i < 4; i ++)
            for(j = 0; j < 4; j ++)
                unit.a[i + 4][j] = Akb.a[i][j];
    }
    for(i = 0; i < 4; i ++)
        mat.a[i][i] = mat.a[i][i + 4] = 1;
    Ak = getA(K);
    for(i = 0; i < 4; i ++)
        for(j = 0; j < 4; j ++)
            mat.a[i + 4][j + 4] = Ak.a[i][j];
    unit = powmod(unit, mat, N - 1);
    printf("%d\n", unit.a[0][0]);
}
int main()
{
    while(scanf("%d%d%d%d", &K, &B, &N, &M) == 4)
    {
        solve();
    }
    return 0;
}