FZU2018级算法第一次作业 1.1fibonacci （矩阵快速幂）

题目

　　Winder最近在学习fibonacci 数列的相关知识。我们都知道fibonacci数列的递推公式是F（n）=F（n-1）+F（n-2）（n>=2 且n 为整数）。 Winder想知道的是当我们将这个递推式改为F（n）=AF（n-1）+BF（n-2）（n>=2且n为整数）时我们得到的是怎样的数列。但是，Winder很懒，所以只能由你来帮他来完成这件事。 注意，这里我们依然令F（0）=F（1）=1。

★数据输入

　　输入第一行三个正整数N，A 和B（N<=10；1<=A、B<=100 且均为整数）。 接下来有N 行，每行一个自然数n（n<=100000000）。

★数据输出

　　输出一行一个整数F（n），由于结果可能会很大，Winder要求输出结果对2013取模。

输入示例	输出示例
5 4 5 2 4 8 16 32	9 209 1377 182 9

题解：

　　一道很经典的矩阵快速幂裸题。

　　首先讲解快速幂，当我们需要求$a^{b}$对mod取模时，可以将b转化为2进制，就可以将b转换为若干个二次幂之和。例如，我们在计算2的12次方时，12的二进制为1100，1100中的两个1的位权分别为4和8，因此$2^{12}$次方便可以转换为$2^{4}*2^{8}$。由此，我们先将答案ans赋值为1，只需要计算a的1，2，4，8.....次方，然后看一下b在该位的数值是否为1即可，如果为1，将ans乘上即可。

快速幂的代码如下

ll fastpow(ll a,ll b)
{
    ll ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}

　　接下来回到本题，由于n的数字较大，加上取模速度较慢，本题递归递推会超时。因此需要寻找复杂度小于O(n)的算法。

　　由于

$
left(
egin{matrix}
a & b \
1 & 0 
end{matrix}
ight)*
inom{f(n)}{f(n-1)}=inom{f(n+1)}{f(n)}
$我们可以构造矩阵，可以得到

$$
left(
egin{matrix}
a & b \
1 & 0
end{matrix}
ight) ag{2}^{n-1}*
inom{f(1)}{f(0)}=inom{f(n)}{f(n-1)}
$$

由此，我们只需要计算矩阵

$$
left(
egin{matrix}
a & b \
1 & 0
end{matrix}
ight)
$$

的n-1次方即可，对于这个矩阵的n-1次方，使用快速幂求出所求矩阵，便可以在$log n$的时间内计算出f(n)的值。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<string>
#include<vector>
#include<queue>
#include<set>
using namespace std;
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
typedef pair <int,int> pii;
#define rep(i,x,y) for(int i=x;i<y;i++)
#define rept(i,x,y) for(int i=x;i<=y;i++)
#define per(i,x,y) for(int i=x;i>=y;i--)
#define pb push_back
#define fi first
#define se second
#define mes(a,b) memset(a,b,sizeof a)
const int inf=0x3f3f3f3f;
const int mod=2013;
class matrix
{
    public:
        int arrcy[6][6];//arrcy为矩阵，下表从0开始
        int row,column;//row为矩阵的行，column为矩阵的列
        friend matrix operator *(matrix s1,matrix s2)
        {
            int i,j;
            matrix s3;
            for (i=0;i<s1.row;i++)
            {
                for (j=0;j<s2.column;j++)
                {
                    for (int k=0;k<s1.column;k++)
                    {
                        s3.arrcy[i][j]+=s1.arrcy[i][k]*s2.arrcy[k][j];
                        s3.arrcy[i][j]%=mod;
                    }
                }
            }
            s3.row=s1.row;
            s3.column=s2.column;
            return s3;
        }
};
matrix quick_pow(matrix s1,long long n)//矩阵快速幂函数，s1为矩阵，n为幂次
{
    matrix mul=s1,ans;
//将ans构造为单位矩阵
    ans.row=ans.column=s1.row;
    memset(ans.arrcy,0,sizeof ans.arrcy);
    for(int i=0;i<ans.row;i++)
        ans.arrcy[i][i]=1;
    while(n)
    {
        if(n&1) ans=ans*mul;
        mul=mul*mul;
        n/=2;
    }
    return ans;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    int n,a,b;
    cin>>n>>a>>b;
    matrix mul;
    mul.row=mul.column=2;
    mul.arrcy[0][0]=a;
    mul.arrcy[0][1]=b;
    mul.arrcy[1][0]=1;
    mul.arrcy[1][1]=0;
    matrix r;
    r.row=2;
    r.column=1;
    r.arrcy[0][0]=r.arrcy[1][0]=1;
    rep(i,0,n)
    {
        int x;
        cin>>x;
        if(!x)//当x=1时，x-1<0无法使用快速幂，答案为0，特判即可
        {
            cout<<1<<endl;
            continue;
        }
        matrix mm=quick_pow(mul,x-1);
        mm=mm*r;
        cout<<mm.arrcy[0][0]<<endl;
    }
    return 0;
}