nyoj1000&&hdu4549 M斐波那契数列 递推公式+矩阵

点击hdu4549

这个题目是矩阵求斐波拉契数列的进阶版，递推关系改成了F(n-1)*F(n-2)=F(n)，所以在用原来的矩阵实现是不可以的，所以这时候就要想别的方法。

既然这个题目是递推公式的题目，所以不妨写出几项来观察规律：

F(n)=F(n-1)*F(n-2)      F(n-1)=F(n-2)*F(n-3)        F(n-2)=F(n-3)*F(n-4)       F(n-3)=F(n-4)*F(n-5)   F(n-4)=F(n-5)*F(n-6)   F(n-5)=F(n-6)*F(n-7)

所以

F(n)= F(n-2)*F(n-3) *F(n-2)   

=>  F(n)=F(n-2)^2*F(n-3)  

=>  F(n)=F(n-3)^3*F(n-4) ^2

=> F(n)=F(n-4)^5*F(n-5) ^3

=> F(n)=F(n-5)^8*F(n-6) ^5 =>.......=>F(n)=F(1)^f(n)*F(0)^f(n-1)

观察上面可得：

F(n)可以成F(n)=F(a)^b*F(a-1)^ c的通式，而且b,c满足斐波拉契数列的规律，所以这个题目需要借助矩阵求出斐波拉契数列的f(n)然后在通过快速幂分别求得F(1)^f(n)与F(0)^f(n-1)

最后相乘。

下面是矩阵求斐波拉契数列的方法

这个题目还需要注意的的是在利用快速幂求F(1)^f(n)与F(0)^f(n-1)的时候，如果对1000000007去取余，那就算用long long任然会溢出，所以这里需要用到数论中的公式结论

费小马定理：

当gcd(A,M)==1的时候 (gcd指最大公约数)
A^X = A^( X mod Eular(M) ) ( mod M ) .这样我们只需对1000000006求余

a^f(n-1) = a^(f(n-1)%1000000006) (mod 1000000007)

b^f(n) = b^(f(n)%1000000006) (mod 1000000007)

这2个直接快速幂就行了

下面是代码：

#include<iostream>
#include<stdio.h>
#include<algorithm>
#include<cstring>
using namespace std;
#define maxn 3
#define me(x) memset(x,0,sizeof(x))
#define ll long long
#define mod1 1000000007
#define mod2 1000000006

struct mat//用结构体封装矩阵
{
    ll m[maxn][maxn];
}unit,temp1,temp2,mm;

void print(mat a)//自己写的打印函数
{
    int i,j;
    for(i=0;i<maxn;i++)
    {
        for(j=0;j<maxn;j++)
            printf("%d ",a.m[i][j]);
        puts("");
    }
    puts("");
    return ;
}
void init()//初始化
{
    me(unit.m);
    me(temp1.m);
    me(temp2.m);
    me(mm.m);

    for(int i=0;i<3;i++)
        unit.m[i][i]=1;

    temp1.m[0][0]=1;temp1.m[0][1]=1;
    temp1.m[1][0]=1;temp1.m[2][1]=1;

    temp2.m[0][0]=1;temp2.m[0][1]=1;
    temp2.m[1][0]=1;temp2.m[2][1]=1;


    return ;
}

mat operator*(mat a,mat b)//对*运算符重载为矩阵相乘
{
    mat p;
    for(int i=0;i<maxn;i++)
    {
        for(int j=0;j<maxn;j++)
        {
            p.m[i][j]=0;
            for(int w=0;w<maxn;w++)
            {
                p.m[i][j]+=a.m[i][w]*b.m[w][j]%mod2;
                p.m[i][j]%=mod2;
            }
        }
    }
    return p;
}
ll pow2(long long a,int b)//¿ìËÙÃÝ Î»ÔËËã°æ
{
    long long r=1,base=a;
    while(b)
    {
    if(b&1) r=(r*base)%mod1;
    base=(base*base)%mod1;
    b>>=1;
    }
    return r%mod1;
}
ll pow(mat N,int k,mat p)//矩阵快速幂
{
    mat ans;
    ll sum=0;
    ans=unit;

    while(k)
    {
        if(k&1)
        {
            ans=ans*N;
        }
        N=N*N;
        k=k>>1;
    }
    return ans.m[0][0]%mod2;//因为ans[0][0]就是f(n);
}

int main()
{
    int a,b,n;
    long long cnt;
    while(scanf("%d %d %d",&a,&b,&n)!=EOF)
    {
        if(n==0||n==1||n==2)
        {
            if(n==0) printf("%d\n",a);
            else if(n==1) printf("%d\n",b);
            else  printf("%d\n",(a*b)%mod1);
            continue;
        }
         init();
         mm.m[0][0]=a*b, mm.m[1][0]=b, mm.m[2][0]=a;
            /*print(mm);
            print(temp1),print(temp2);
            cout<<n-2<<" "<<n-3<<endl;*/

         ll a1=pow(temp1,n-1,mm),a2=pow(temp2,n-2,mm);///其中a,b为M斐波拉契数列中第0,1项，所以第n-1项相当于n项

         // printf("%lld %lld\n",a1,a2);

         a1=pow2(b,a1),a2=pow2(a,a2);
         printf("%lld\n",(a1*a2)%mod1);
         
    }
    return 0;
}