hdu1588 Gauss Fibonacci 矩阵快速幂

$F(x)$是斐波拉契数列
给定一次函数$g(x)=ax+b$
求 $\underset {i=0} {\overset n {\Large \Sigma}} F(g(i))$
抄来的latex的写法，哦哟这个公式看起来逼格就是高

根据矩阵乘法的思想，我们可以把$F(n)=F(n-1)+F(n-2)$ 作成变换矩阵 $R$ ，连乘a次后可以得到一个变换矩阵$R^a$，那么可以用它做出如下变换：

$${R^a}* \begin {pmatrix} F(n)\\ F(n-1) \end {pmatrix} {=} \begin {pmatrix} F(n+a)\\ F(n+a-1) \end {pmatrix}$$

然而我们要求的是若干斐波拉契数的和
所幸的是所求项是有规律的，是从b开始每隔a选中的一个数

上面的矩阵乘法没有保存数，那我们多开个坑保存不就行了
原向量：

$$\begin {pmatrix} F(n)\\ F(n-1) \end {pmatrix}$$
修改后：
$$\begin {pmatrix} F(n)\\ F(n-1)\\ Sum \end {pmatrix} (1)$$
然后我们修改变换矩阵，使得它能够完成把$F(n-1)$累加在 $Sum$中（为什么不累加$F(n)$？因为感觉会出来一堆边界判断，太傻）
原斐波拉契数列变换矩阵

修改后（左上角四个$?$组成$R^a$)

$$\begin {pmatrix} ?&?&0\\ ?&?&0\\ 0&1&1\\ \end {pmatrix} (2)$$

现在，拿修改后的矩阵去左乘(1)，即可得到

$$\begin {pmatrix} F(n+a)\\ F(n+a-1)\\ Sum+F(n-1) \end {pmatrix} (3)$$

这就是 我要的 滑板鞋，哟 哟 哟

回归到题目，给定一次函数$g(x)=ax+b$，先建立2*1的向量即F(0)和F(1)，和2*2的小矩阵R，即普通斐波拉契变换的矩阵

$$初始向量： \begin{pmatrix} F(1)\\ F(0)\\ \end{pmatrix} 斐波拉契变换矩阵R： \begin{pmatrix} 1&1\\ 1&0\\ \end{pmatrix}$$
快速幂求出$R^{b}$，左乘初始向量得到如下向量
$$\begin{pmatrix} F(b+1)\\ F(b)\\ 0 \end{pmatrix} (4)$$
再快速幂求出$R^a$，扩展成3*3，修改为$(2)$ 的形式，记作$R_3$，
求出${R_3}^{n}$ 然后左乘向量(4)即可
得到向量的Sum位置放着$F(b)+F(b+a)+F(b+2*a)+……F(b+n*a)$，就是最后的答案啦

直接把上题的模板套来用了，我好强啊
main中的内容和步骤是一模一样的

注意long long

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<map>
#include<vector>
#include<cmath>
using namespace std;


#define ll long long
//Template mat
//Interface here
#define MULMOD
//Use mat(int row,int cow) to define a new EMPTY mat
//Use define MULMOD and set MOD to ENABLE MOD
#define MAXROWS 4
#define MAXCOLS 4
#ifdef MULMOD
int MOD;
#endif
//Interface End

struct Matrix
{
    int Rows,Cols;
    ll data[MAXROWS][MAXCOLS];
    void clear()
    {
        memset(data,0,sizeof(data));
    }
    Matrix (int n,int m)
    :Rows(n),Cols(m)
    {
        clear();
    }
    Matrix (int n)
    :Rows(n),Cols(n)
    {
        clear();
        for (int i=0;i<n;i++)
            data[i][i]=1;
    }
    ll* operator[](const int n)
    {
        return data[n];
    }
    Matrix operator* (const Matrix& ano) const
    {
        Matrix result(Rows,ano.Cols);
        for (int i=0;i<Rows;i++)
            for (int j=0;j<ano.Cols;j++)
                for (int k=0;k<Cols;k++)
                {
                    #ifdef MULMOD
                    result[i][j] += data[i][k] * ano.data[k][j] % MOD;
                    result[i][j]%=MOD;
                    #else
                    result[i][j] += data[i][k] * ano.data[k][j];
                    #endif
                }
        return result;
    }
};
Matrix QuickMatrixPow(Matrix to ,int k)
{
    Matrix ans(to.Rows);
    while (k)
    {
        if (k&1)
            ans=ans*to;
        k>>=1;
        to = to*to;
    }
    return ans;
}
//Template mat end


int main()
{
    cin.sync_with_stdio(false);
    int a,b,n,m;
    while (cin>>a>>b>>n>>m)
    {
        MOD=m;
        Matrix V0(2,1);
        V0[0][0]=1;
        V0[1][0]=0;
        Matrix R(2,2);
        R[0][0]=1;
        R[0][1]=1;
        R[1][0]=1;
        R[1][1]=0;
        Matrix Vb(QuickMatrixPow(R,b)*V0);
        Vb.Rows=3;      //(4)
        Vb[2][0]=0;     //Sum
        Matrix Ra(QuickMatrixPow(R,a));
        Ra.Rows=3;
        Ra.Cols=3;
        Ra[2][1]=Ra[2][2]=1;
        Ra[0][2]=Ra[1][2]=Ra[2][0]=0;
        cout<< (QuickMatrixPow(Ra,n) * Vb )[2][0] <<endl;
    }
}