矩阵快速幂（二）

浅谈矩阵快速幂，希望能够帮助到你！

现在给出一个实际的问题，要求2的十次方，很简单，直接for循环就行了。

可是现在如果说让求2的10000000000次方，又该怎么办呢？

这个时候我们可以考虑使用快速幂。

#include<algorithm>
#include<iostream>
#include<string.h>
#include<cstdio>
#define  P pair<long long,long long>
#define ll long long
using namespace std;
ll x,y;
ll pow(ll x,ll y)
{
    ll ans=1;
    while(y)
    {
        if(y%2==1)
        {
            ans*=x;
        }
        x*=x;
        y>>=1;
    }
    return ans;
}

int main()
{
        cin>>x>>y;
        printf("%lld\n",pow(x,y));
    return 0;
}

矩阵快速幂分为两个部分：

1：奇数

2：偶数

 

从100直接到50，就相当于二分，所以能够大大降低时间复杂度。

 

现在把问题从快速幂转移到矩阵快速幂。

 

第一个问题是求斐波那契数：

在Fibonacci整数序列中，对于n≥2，F0 = 0，F1 = 1，Fn = Fn  -  1 + Fn  -  2.例如，Fibonacci序列的前十项是：
0,1,1,2,3,5,8,13,21,34 ......

给定整数n，您的目标是计算Fn的最后4位数。
输入
输入测试文件将包含多个测试用例。每个测试用例由一行包含n（其中0≤n≤1,000,000,000）组成。文件结尾由包含数字-1的单行表示。
输出
对于每个测试用例，打印Fn的最后四位数字。如果Fnare的最后四位数全为零，则打印'0';否则，省略任何前导零（即，打印Fn mod 10000）。
样本输入
0
9
999999999
-1
样本输出
0
34
626

要求的是斐波那契数，我们知道对于斐波那契数的第n位（n>=3）满足f(n)=f(n-1)+f(n-2)。

可以推出A是：

也就是说，现在我们要求F(N)，只需要对A进行N-1次方就行了。

具体的实现如下：

#include<algorithm>
#include<iostream>
#include<string.h>
#include<cstdio>
#include<cstring>
#define  P pair<long long,long long>
#define ms(a,b) memset(a,b,sizeof(a))
#define ll long long
using namespace std;
const int mod=1e4;
//ll x,y;
//ll pow(ll x,ll y)
//{
//    ll ans=1;
//    while(y)
//    {
//        if(y%2==1)
//        {
//            ans*=x;
//        }
//        x*=x;
//        y>>=1;
//    }
//    return ans;
//}
struct xiao
{
    ll f[4][4];
} ans,mapp;

void init(xiao &ans)
{
    for(int i=1; i<=2; ++i)
    {
        for(int j=1; j<=2; ++j)
        {
            if(i==j)
                ans.f[i][j]=1;
            else
                ans.f[i][j]=0;
        }
    }

}

void power(xiao &x,xiao &y,xiao &z)
{
    //z初始化，这里是存储相乘后矩阵信息的矩阵
    ms(z.f,0);
    for(int i=1; i<=2; ++i)
    {
        for(int j=1; j<=2; ++j)
        {
            if(x.f[i][j])
            {
                for(int k=1; k<=2; ++k)
                {
                    z.f[i][k]+=x.f[i][j]*y.f[j][k];
                    z.f[i][k]%=mod;
                }
            }
        }
    }
    /*
    注意，这里把所有的信息都存储到z里面了，现在的情况只是说我有100毫升水在C杯子里面，
    我往A中倒了50ml，现在还有50ml在c杯子里面，接下来就是把剩下的50ml水倒到b杯子里面。
    实现这一步就是简单的赋值。具体实现在下面.

    */
}

ll judex(ll x)
{
    init(ans);
    /*要求第n位，对n-1位进行快速幂运算。
      类比前面提到的快速幂的知识，我们首先定义一个单位矩阵，相当于前面的那个ans变量。
    */
    xiao temp=mapp,t;
    x--;
    while(x)
    {
        if(x%2==1)
        {
            //类比前面的 ans*=2;
            power(ans,temp,t);
            ans=t;
        }
        //这里类比把2的100次方变成4的50次方。
        power(temp,temp,t);
        temp=t;
        x/=2;
    }
    return ans.f[1][2]+ans.f[1][1];
}

int main()
{
    mapp.f[1][1]=0;
    mapp.f[1][2]=1;
    mapp.f[2][1]=1;
    mapp.f[2][2]=1;

    ll n;//要求的是哪一位的
    while(cin>>n)
    {
        if(n==-1)
            break;
        if(n==1)
        {
            printf("1\n");
            continue;
        }
        else if(n==2 )
        {
            printf("1\n");
            continue;
        }
        else if(n==0)
        {
            printf("0\n");
            continue;
        };
        printf("%lld\n",judex(n)%mod);
    }
    return 0;
}

 

如果这道题解决了的话，接下来再看一道稍微复杂一点的题：

Recursive sequence

 

农夫约翰喜欢和他的N奶牛一起玩数学游戏。最近，它们被递归序列所吸引。在每个回合中，奶牛会排成一排，而约翰在黑板上写下两个正数a和b。然后，奶牛会逐一说出他们的身份证号码。第一头牛说第一个数字a，第二头说第二个数字b。在那之后，第i个牛说出第（i-2）个数，第（i-1）个数和i4的两倍之和。现在，您需要编写一个程序来计算第N头奶牛的数量，以检查John的奶牛是否可以做到正确。

输入

第一行输入包含一个整数t，即测试用例的数量。 t测试案例如下。 每个情况仅包含一行具有三个数字N，a和b，其中N，a，b <231，如上所述。

输出

对于每个测试用例，输出第N个牛的数量。这个数字可能非常大，所以你需要以模数2147493647输出它。

样本输入

2

3  1  2

4  1  10

样例输出

85

369

 

接下来，我们一起来一步步分析一下这道题：

首先，我们知道的是F(n)=2*f(n-2)+f(n-1)+n^4.

那么F(n+1)=2*f(n-1)+f(n)+（n+1）^4.

现在的思想还是类比的思想，我们在前面的斐波那契数的计算的时候知道，使用矩阵实现求得f(n)，那么我们也可以使用矩阵实现从n^4到(n+1)^4的实现！

接下来请看下面的图：

也就是说，既然可以拆分成n^4,n^3,n^2,n^1,1这些项，那这些不就可以作为矩阵的内容吗？

对，那就可以实现了！接下来就是求A了。

因为矩阵里面有7个元素，所以我们构造的矩阵是7*7的矩阵。

举个例子，对于(n+1)^4来说，拆开就是

n^4+4*n^3+6*n^2+4*n+1.

按照这样的方法可以求出A

A=

构造算是矩阵快速幂的难点了，构造出来以后就好做了，接下来的一个难点是输出。

我们知道，输入的三个数分别为k,a,b，现在k如果是1的话，输出的是a，如果k是2的话，输出的是b，如果k是三的话，输出的就是2*a+b+1,这样就开始使用矩阵乘了。

    return ((ans.f[1][0])*a+(ans.f[1][1])*b+(ans.f[1][2])*16+(ans.f[1][3])*8+(ans.f[1][4])*4+(ans.f[1][5])*2+(ans.f[1][6]));

附上完整的代码：

#include<algorithm>
#include<iostream>
#include<string.h>
#include<cstdio>
#include<cstring>
#define  P pair<long long,long long>
#define ms(a,b) memset(a,b,sizeof(a))
#define ll long long
using namespace std;
const ll mod=2147493647;
ll N,a,b,k;

struct xiao
{
    ll f[7][7];
} ans,mapp;

void init(xiao &ans)
{
    for(int i=0; i<7; ++i)
    {
        for(int j=0; j<7; ++j)
        {
            if(i==j)
                ans.f[i][j]=1;
            else
                ans.f[i][j]=0;
        }
    }
}

void power(xiao &x,xiao &y,xiao &z)
{
    ms(z.f,0);
    for(int i=0; i<7; ++i)
    {
        for(int j=0; j<7; ++j)
        {
            if(x.f[i][j])
            {
                for(int k=0; k<7; ++k)
                {
                    z.f[i][k]+=(x.f[i][j]*y.f[j][k]);
                    z.f[i][k]%=mod;
                }
            }
        }
    }
}

ll solve(ll a,ll b,ll k)
{
    init(ans);
    ll number=k-2;
    xiao temp=mapp,t;
    while(number)
    {
        if(number%2==1)
        {
            power(ans,temp,t);
            ans=t;
        }
        power(temp,temp,t);
        temp=t;
        number/=2;
    }
    return ((ans.f[1][0])*a+(ans.f[1][1])*b+(ans.f[1][2])*16+(ans.f[1][3])*8+(ans.f[1][4])*4+(ans.f[1][5])*2+(ans.f[1][6]));
}


int main()
{
    mapp=
    {
        0,1,0,0,0,0,0,
        2,1,1,4,6,4,1,
        0,0,1,4,6,4,1,
        0,0,0,1,3,3,1,
        0,0,0,0,1,2,1,
        0,0,0,0,0,1,1,
        0,0,0,0,0,0,1
    };
    for(cin>>N; N--;)
    {
        cin>>k>>a>>b;
        if(k==1)
        {
            printf("%lld\n",a);
        }
        else if(k==2)
        {
            printf("%lld\n",b);
        }
        else
        printf("%lld\n",solve(a,b,k)%mod);
    }
    return 0;
}

 

最后针对矩阵的乘法进行补充：
 

最后得到的矩阵的形状是：

那么我们按照遍历一遍的方法去做的话，我们不能说每一个坐标点我们只计算一次，但我们可以肯定的是每一行我们只走一次，那么我们在走过这一行的时候就要把所有的数值都给计算出来。

对于上面的思想，我们用到的方法是矩阵乘法里面的把每一行都给提取出来进行计算（前面矩阵中的行）。