POJ 3040 浅谈矩阵快速幂优化LogN斐波拉契函数求解

世界真的很大
讲道理矩阵这个东西还完全不怎么会怎么办？
赶紧复习一波再找一道水题刷一刷
矩阵优化线性函数求解曾经也听说过，但今天不知怎么奇困无比然后公式推起来巨慢
但其实还是比较简单

看题先：

description：

求斐波拉契数列第n项

input：

The input test file will contain multiple test cases. Each test case consists of a single line containing n (where 0 ≤ n ≤ 1,000,000,000). The end-of-file is denoted by a single line containing the number −1.

output：

For each test case, print the last four digits of Fn. If the last four digits of Fn are all zeros, print ‘0’; otherwise, omit any leading zeros (i.e., print Fn mod 10000).

讲道理斐波拉契函数本来就可以O（n）求出来，看上去已经很优了
但是这道题的数据范围奇大无比，1e9
基于斐波拉契是一个线性函数，而且每一项的值只与前两项有关
而且是和的关系
只要碰上这样的函数，可以考虑由矩阵快速幂优化求解

由于答案与前两项有关，那么我们很自然想到搞一个只有两项的矩阵。一个1*2的矩阵a，a（1,1）是F（i），a（2,1）是F（i-1 )
这样就很轻松得到转移矩阵e
（ 1， 1 ）
（ 1， 0 ）
求F（n）就等于F（1）再乘以e^(n-1)，这一步用快速幂搞定

切记注意矩阵乘法的运算顺序，这是由于其没有交换律的缘故

完整代码：

#include<stdio.h>
#include<cstring>
using namespace std;
typedef long long dnt;

const int mod=10000;

struct Matrix
{
    int x,y;
    dnt a[5][5];
    void reset()
    {
        memset(a,0,sizeof(a));
    }
}e,F;

Matrix operator*(const Matrix &aa,const Matrix &bb)
{
    Matrix rt;
    rt.reset();
    rt.x=aa.x,rt.y=bb.y;
    for(int i=1;i<=rt.x;i++)
        for(int j=1;j<=rt.y;j++)
            for(int k=1;k<=aa.y;k++)
                rt.a[i][j]=(rt.a[i][j]+aa.a[i][k]*bb.a[k][j])%mod;
    return rt;
}

Matrix quickmub(Matrix a,int b)
{
    Matrix rt;
    rt.x=rt.y=a.x;
    rt.reset();
    for(int i=1;i<=rt.x;i++) rt.a[i][i]=1;
    while(b)
    {
        if(b&1) rt=rt*a;
        a=a*a;
        b>>=1;
    }
    return rt;
}

int main()
{
    F.x=2,F.y=1,F.reset();
    F.a[1][1]=1,F.a[2][1]=0;
    e.x=2,e.y=2,e.reset();
    e.a[1][1]=1,e.a[1][2]=1,e.a[2][1]=1,e.a[2][2]=0;
    while(1)
    {
        int n;
        scanf("%d",&n);
        if(n==-1) break ;
        Matrix A;
        A.reset();
        if(n>0) A=quickmub(e,n-1)*F;
        printf("%d\n",A.a[1][1]);
    }
    return 0;
}
/*
Whoso pulleth out this sword from this stone and anvil is duly born King of all England
*/

嗯，就是这样