M斐波那契数列

题目描述：
M斐波那契数列F[n]是一种整数数列，它的定义如下：
F[0] = a
F[1] = b
F[n] = F[n-1] * F[n-2] ( n > 1 )
现在给出a, b, n，你能求出F[n]的值吗？（ 0 <= a, b, n <= 10^9 ）
由于F[n]可能很大，你只需输出F[n]对1000000007取模后的值即可

题解：
列出F[n]的前几项会发现，F[n] = a^f(n-2) * b^f(n-1)，f(n)为斐波拉契数列，用矩阵快速幂求解
费马小定理：a^k % p = a^(k % (p-1)) 其中，p 为质数
a^f(n-2) * b^f(n-1) % p = a^( f(n-2) % (p-1) ) * b^( f(n-1) % (p-1) ) % p

代码

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <iostream>
#include <algorithm>
#include <string>
#include <set>
#include <map>
#include <queue>
#include <stack>
#include <vector>
using namespace std;

const int mod = 1e9+7;
int n,x,y;

struct Matrix
{
    int n,m,d[5][5];
    Matrix (int a,int b)
    {
        n = a, m = b;
        memset(d,0,sizeof(d));
    }
    void Copy(int *tmp)
    {
        for (int i=0;i<n;i++)
            for (int j=0;j<m;j++)
            {
                d[i][j] = (*tmp);
                tmp++;
            }
    }
    friend Matrix operator * (const Matrix &a,const Matrix &b)
    {
        Matrix c(a.n,b.m);
        for (int i=0;i<a.n;i++)
            for (int j=0;j<b.m;j++)
            {
                long long tmp = 0;
                for (int k=0;k<a.m;k++)
                {
                    tmp += (long long)a.d[i][k] * b.d[k][j] % (mod-1);
                    tmp %= (mod-1);
                }
                c.d[i][j] = tmp % (mod-1);
            }
        return c;
    }
};

long long qpow(int x,int k)
{
    long long ans = 1;
    x %= mod;
    while (k)
    {
        if (k & 1) ans = (ans*x) % mod;
        x = (long long)x*x % mod;
        k >>= 1;
    }
    return ans;
}

long long solve(int x,int y,int n)
{
    if (n == -1) return x;
    if (n == 0) return y;
    Matrix A(1,2),B(2,2);
    int a[2] = {1,0};
    int b[4] = {1,1,1,0};
    A.Copy(a);
    B.Copy(b);
    while (n)
    {
        if (n & 1) A = A*B;
        B = B*B;
        n >>= 1;
    }
    long long f=0,g=0;
    f=qpow(x,A.d[0][1]);
    g=qpow(y,A.d[0][0]);
    return (f * g % mod);
}

int main()
{
    while (scanf("%d %d %d",&x,&y,&n)!=EOF)
    {
        printf("%I64d\n",solve(x,y,n-1));
    }
    return 0;
}