[POJ 1588] Gauss Fibonacci （矩阵快速幂）

链接

HDU 1588

---

题意

求若干项fibonacci数列的sum(fibonacci[j])，j = k*i + b。
fibonacci数列的等差数列项求和。

---

题解

矩阵快速幂的一个经典例题就是快速求出某项fibonacci数。
还有一个经典问题，可解决A + A^2 + … + A^n，矩阵快速幂求和。
事实上这道题目是两者的综合，fibon的第k项等价为某矩阵F的k次幂F^K（中的某项），第2k项等价为F^2K，第nK项等价为F^nK。。。
设F^K为A，求A + A^2 + … + A^n即可，可以看到其实是有数学推理在里面的。

细节要注意乘法溢出，0边缘的各数试了几次都通过了，提交就AC了。

---

代码

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <iostream>
using namespace std;
typedef long long lint;
int mod;
struct matrix
{
    int a[4][4], n;
    void clear() { memset(a, 0, sizeof(a)); }
    matrix(int k, int type)
    {
        n = k, clear();
        if(type) for(int i = 0; i < n; i++)
            a[i][i] = 1;
    }
    matrix() { n = 2, clear(); a[0][0] = a[0][1] = a[1][0] = 1, a[1][1] = 0; }
    matrix operator* (const matrix& b) const
    {
        matrix o = matrix(n, 0);
        for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++)
        for(int k = 0; k < n; k++)
        {
            o.a[i][j] += ((lint)a[i][k] * b.a[k][j]) % mod;
            o.a[i][j] %= mod;
        }
        return o;
    }
    friend matrix operator^ (matrix tmp, int k)
    {
        matrix o = matrix(tmp.n, 1);
        while(k)
        {
            if(k & 1) o = o * tmp;
            tmp = tmp * tmp;
            k >>= 1;
        }
        return o;
    }
};
matrix getKthsum(const matrix& b, int k)
{
    matrix tmp = matrix(b.n << 1, 0);
    for(int i = 0; i < b.n; i++)
    for(int j = 0; j < b.n; j++)
    tmp.a[i][j] = b.a[i][j];
    for(int i = 0; i < b.n; i++)
    {
        tmp.a[i][b.n + i] = 1;
        tmp.a[b.n + i][b.n + i] = 1;
    }
    tmp = tmp^(k+1);
    matrix o = matrix(b.n, 0);
    for(int i = 0; i < b.n; i++)
    for(int j = 0; j < b.n; j++)
    o.a[i][j] = tmp.a[i][b.n + j];
    for(int i = 0; i < b.n; i++)
        o.a[i][i] = (o.a[i][i] + mod - 1) % mod;

    return o;
}
int main()
{
    int k, b, n, M;
    while(cin >> k >> b >> n >> M)
    {
        mod = M;
        matrix fmat = matrix();
        matrix dmat = fmat ^ b;
        matrix base = fmat ^ k;
        matrix ans = getKthsum(base, n - 1);
        ans = ans * dmat;
        cout << (ans.a[0][1] + dmat.a[0][1]) % mod << endl;
    }
    return 0;
}
//1000000000