HDU 1588 矩阵快速幂 嵌套矩阵

这个题目搞了我差不多一个下午，之前自己推出一个公式，即 f[n+k]=k*f[n]+f[n-1]结果发现根本不能用，无法降低复杂度。

后来又个博客的做法相当叼，就按他的做法来了

即 最终求得是 S（n）=f[b]+f[b+k]+f[b+2*k]....f[b+n*k] (原题的意思好像是不用加到第n项，但实测确实要加到该项)

然后我们令 A={1，1}（标准的斐波那契矩阵）

                    {1,0}发现 f[b]=A^b,f[b+k]=A^(b+k),....f[b+nk]=A^(b+nk);

提取公共因子 A^b.S(n)=A^b*(E+A^K+A^K^2....A^K^n)

再令K=A^K  (K和E都是矩阵，E为单位矩阵)。于是S（n）=A^b*(E+K+k^2+K^3+...K^n);

这样我们的目的大概出来了，就是要尽可能短的算出 上式括号中的内容(A^b可以用经典斐波那契矩阵很快算出来)

于是构造一个嵌套矩阵 

令 SK={K，E}，（K,E,0均为矩阵），SK^N={K^N ,  E+K+K^2+K^3+...K^N}

           {0，E}                                         {0    ,  E       }

这样只要求出SK^N,取其第一行的第二项 再与 A^b相乘，即可得出最终结果

嵌套矩阵跟普通矩阵差不多，调用写好的矩阵乘法即可实现嵌套矩阵的乘法和乘方

#include <iostream>
#include <cstdio>
#include <cstring>
#define ll __int64
using namespace std;
ll k,b,n,M;
struct Mat
{
    ll mat[3][3];
} Fb,E,zero,A;
struct Smart
{
    Mat d[3][3];
} SK,SE;
void test(Mat t)
{
    for (int i=0;i<2;i++)
    {
        for (int j=0;j<2;j++)
            cout<<t.mat[i][j]<<" ";
        cout<<endl;
    }
}
void test2(Smart tmp)
{
    for (int i=0;i<2;i++)
    {
        for (int j=0;j<2;j++)
            test(tmp.d[i][j]);
    }
}
Mat operator +(Mat a,Mat b)
{
    Mat 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];
            c.mat[i][j]%=M;
        }
    }
    return c;
}
Mat operator *(Mat a,Mat b)
{
    Mat c;
    memset(c.mat,0,sizeof c.mat);
    for (int i=0;i<2;i++)
    {
        for (int j=0;j<2;j++)
        {
            for (int k=0;k<2;k++)
            {
                c.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
                if (c.mat[i][j]>=M)
                    c.mat[i][j]%=M;
            }
        }
    }
    return c;
}
Mat operator ^(Mat a,ll x)
{
    //if (x==0) return a;
   // cout<<x<<" the x"<<endl;
    Mat c=E;
    for (;x;x>>=1)
    {
        if (x&1)
            c=c*a;
        a=a*a;
    }
    return c;
}
Smart operator *(Smart a,Smart b) //大矩阵的乘法 调用普通矩阵的乘法和加法公式即可
{
    Smart c;
    Mat tmp;
    for (int i=0;i<2;i++)
    {
        for (int j=0;j<2;j++)
        {
            c.d[i][j]=zero;
            for (int k=0;k<2;k++)
            {
                tmp=a.d[i][k]*b.d[k][j];
                c.d[i][j]=c.d[i][j]+tmp;
            }
        }
    }
    return c;
}
Smart operator ^(Smart a,ll x) //大矩阵的乘法 调用普通矩阵的乘法公式即可
{
    Smart c=SE;
   // if (x==0) return a;
    //cout<<x<<" the s x "<<endl;
    for (;x;x>>=1)
    {
        if (x&1)
            c=c*a;
        a=a*a;
      // test2(c);
    }
   // test2(c);
    return c;
}
void init()
{
    memset(zero.mat,0,sizeof zero);
    memset(E.mat,0,sizeof E.mat);
    for (int i=0;i<2;i++)
        E.mat[i][i]=1;
    Mat A;
    memset(A.mat,0,sizeof A.mat);

    SE.d[0][0]=SE.d[1][1]=E;
    SE.d[0][1]=SE.d[1][0]=zero;

    A.mat[0][0]=A.mat[0][1]=A.mat[1][0]=1;
    Fb=A^b;
    //test(Fb);
    SK.d[0][0]=A^k;
    SK.d[0][1]=SK.d[1][1]=E;
    SK.d[1][0]=zero;
}
int main()
{

    while (scanf("%I64d %I64d %I64d %I64d",&k,&b,&n,&M)!=EOF)
    {
        init();
        Smart s=SK^n;
        Mat tmp=s.d[0][1];
        //test(tmp);
        Mat ans=Fb*tmp;
        printf("%I64d
",ans.mat[1][0]);//因为Fb求出来 真正的 A^b在[1][0]处，但是最终结果为何也是1 0处 我觉得还是有待考究

    }
return 0;
}