矩阵加速——斐波那契数列

来自洛谷P1962(一道看似很水的题)

斐波那契数列的通项公式是 Fn=Fn-1 + Fn-2

在一定的复杂度内可以直接递推,但是如果n太大,那么就容易T,这时候,我们就运用矩阵加速来进行优化,以减少运行时间。


在看矩阵加速之前,我们首先需要了解矩阵快速幂

【模板】 洛谷P3390

 

 首先,我们来讲一下矩阵与矩阵之间的运算。

1.矩阵加法:

假定有两个矩阵A,B;

一般而言,让我们进行矩阵加法的两个矩阵会是一对同型矩阵(行列数分别相等)

那么,得到的矩阵AB的第 i 行,第 j 列的元素,即为A和B中第 i 行,第 j 列的元素的加和

举个栗子:

 那么减法运算就相当于加上一个负数,与之一样。

在加(减)法中,满足:

A+B=B+A【交换律】

A+(B+C)=(A+B)+C 【结合律】

 

2.矩阵数乘

假定我们有一个矩阵A,还有一个常数 λ,那么 λ*A 就是把A中的每一个元素都乘以 λ

再举个例子:

 在矩阵的数乘中,满足如下运算规律:

设 x,y为常数

则:x(yA)= y(xA)【交换律】

     (x+y)*A = xA + yA  【分配律】

     (A+B)*x = xA + xB  【分配律】

3.矩阵乘法(矩阵乘矩阵)

假定我们有两个矩阵A,B。A为n*m的矩阵,B为m*k的矩阵,只有满足A的行数数等于B的列数时,才能够进行矩阵与矩阵之间的乘法运算。运算后,得到的矩阵C将会是一个n*k大小的矩阵

那么得到C矩阵的运算公式为:

 

 又双叒叕举个例子:

 

 

就是说,矩阵AB的第 i 行第 j 列元素都是由A的第 i 行的各个元素与B的第 j 列的各个元素相乘之和

需要注意的是,在矩阵与矩阵的乘法中,有如下运算法则

1.(AB)C=A(BC)【结合律】

2.(A+B)C=AC+BC   【分配律】

这里并没有交换律,因为将A和B的顺序颠倒后,不一定符合矩阵之间的运算定义(比如A:2*3 ; B:3*5),如果满足也不一定能够得到相同的结果。

我们在这道模板题中,因为A只是在做幂运算,并不用加入其他的递推式,我们就可以设定一个单位矩阵 L(单位矩阵的元素非0即1),使得矩阵L的第 i 行,第 i 列的元素为1,其余为0。之后可以依照快速幂的写法完成矩阵快速幂(这里用的重载运算符)

#include<iostream>
#include<cstdio>
#include<cstring>
#include<math.h>
#include<algorithm>
using namespace std;
const long long mod=1000000007;

struct STU //声明结构体
{
    long long m[110][110];
} A,I;

long long n,k;

STU operator * (const STU &a,const STU &b)//重载运算符 '*'
{
    STU x;
    
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            x.m[i][j]=0;
        }
    }
    
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            for(int p=1;p<=n;p++)
            {
                x.m[i][j]+=a.m[i][p]*b.m[p][j]%mod; //矩阵乘法运算
                x.m[i][j]%=mod;//别忘记取mod
            }
        }
    }
    
    return x;
}

int main(void)
{
    cin>>n>>k;
    
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            cin>>A.m[i][j];//输入运算矩阵
        }
        
        I.m[i][i]=1;//定义单位矩阵
    }
    
    while(k>0)//快速幂运算(可以用while循环,也可以用递归)
    {
        if(k%2==1) I=I*A;
        A=A*A;
        k=k>>1;
    } 
    
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            cout<<I.m[i][j]<<" ";//输出每一个元素
        }
        
        cout<<endl;
    }
    
    return 0;
}



既然我们解决了矩阵快速幂的问题,那么再来看一下矩阵加速——【模板】洛谷 P1939(和斐波那契数列的做法差不多)。

 对于这道题来说,直接递归毫无疑问是要T的,那么矩阵加速就是我们的首选算法,那么,问题来了,怎么设置我们的初始矩阵呢?

 在这道题中,我们需要用第n-1个数和第n-3个数来推得第n个数

那么设现有数Fn,Fn-1,Fn-2,Fn-3

那么:Fn=Fn-1*1+Fn-2*0+Fn-3*1 ;Fn-1=Fn-1*1+Fn-2*0+Fn-3*0 ; Fn-2=Fn-1*0+Fn-2*1+Fn-3*0

则矩阵为:

                    

那么这就是我们的初始矩阵。

值得注意的是,在这个矩阵中,第N次运算得到的结果为第N+1项,那么我们就需要输出第二行,第一列的元素即可。

上代码:

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<math.h>
#include<cstring>
using namespace std;

const long long mod=1000000007;

struct STU  //结构体
{
    long long m[110][110];
} A,I;

int T,n;

STU operator * (const STU &a,const STU &b)//重载运算符
{
    STU x;
    
    for(int i=1;i<=3;i++)
    {
        for(int j=1;j<=3;j++)
        {
            x.m[i][j]=0;
        }
    }
    
    for(int i=1;i<=3;i++)
    {
        for(int j=1;j<=3;j++)
        {
            for(int p=1;p<=3;p++)
            {
                x.m[i][j]+=a.m[i][p]*b.m[p][j]%mod;
                x.m[i][j]%=mod;
            }
        }
    }
    
    return x;
}

void clear()
{
    memset(I.m,0,sizeof(I.m));//清空I矩阵
    memset(A.m,0,sizeof(A.m));//清空A矩阵
    for(int i=1;i<=3;i++) I.m[i][i]=1;//还原为初始状态
    A.m[1][1]=A.m[1][3]=A.m[2][1]=A.m[3][2]=1;
} 

int ksm(int k)//快速幂
{
    while(k>0)
    {
        if(k%2!=0) I=I*A;
        A=A*A;
        k=k>>1;
    }
    
    return I.m[2][1];
}

int ys(int k)//判断
{
    if(k<=3) return 1;
    else return ksm(k);
}

int main(void)
{
    cin>>T;
    
    for(int i=1;i<=T;i++)
    {
        cin>>n;
        
        clear();//每一次运算后都需要清空
        
        cout<<ys(n)<<endl;
    }
    
    return 0;
}



正题终于要来了!

 看过了上面的矩阵加速例题,这道题也就变得迎刃而解。

首先,我们先定义初始矩阵,因为它只有前两项为1,那么就可以用一个1*2的矩阵L来代表,即为L【1,1】

之后,我们要定义单位矩阵,由于Fn=Fn-1+Fn-2,那么单位矩阵的第一列就应该定为:【1,1】,使得Fn-1和Fn-2能够相加

之后,为了保留下Fn-1,但Fn-2已经没有用处了,那么第二列就是【1,0】。

这样,我们的单位矩阵就是:

                                                

原式也可以化为:

                        

到这里,我们的矩阵就构造完成了,输出时只要输出第1行第1列的元素就行了,但是,还有一个地方,那就是当我们输入N,且N大于等于3时,我们的矩阵只需要进行N-2运算,因为在本题中第1项和第二项在初始化时已经给出,进行运算时可以直接从求第三项开始。

源代码如下:

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<math.h>
using namespace std;

const long long mod=1000000007;

struct STU//结构体
{
    long long m[110][110];
} A,I;

long long n;

STU operator * (const STU &a,const STU &b)//重载运算符
{
    STU x;
    
    memset(x.m,0,sizeof(x.m));
    
    for(int i=1;i<=2;i++)
    {
        for(int j=1;j<=2;j++)
        {
            for(int k=1;k<=2;k++)
            {
                x.m[i][j]+=a.m[i][k]*b.m[k][j]%mod;
                x.m[i][j]%=mod;
            }
        }
    }
    
    return x;
}
void init()//初始化
{
    memset(I.m,0,sizeof(I.m));
    memset(A.m,0,sizeof(A.m));
    I.m[1][1]=I.m[1][2]=1;
    A.m[1][1]=A.m[1][2]=A.m[2][1]=1;
}
int main(void)
{
    cin>>n;
    
    init();
    
    if(n<=2)//第一项或第二项的情况
    {
        cout<<1<<endl;
        
        return 0;
    }
    else //n>=3的情况
    {
        n=n-2; 
        
        while(n>0) //快速幂
        {
            if(n%2!=0) I=I*A;
            A=A*A;
            n=n>>1;
        }
        
        cout<<I.m[1][1]<<endl;//输出
        
        return 0;
    }
} 

 这样,斐波那契数列的矩阵加速解法就写完啦!

QAQ

原文地址:https://www.cnblogs.com/jd1412/p/12724329.html