hdu3802 矩阵快速幂&特征方程&降幂

转载请注明出处，谢谢http://blog.csdn.net/bigtiao097?viewmode=contents

题意：

给定a、b、p、n,求以下式子的值

$$G(a,b,n,p) = (a^{\frac{p-1}{2}}+1)(b^{\frac{p-1}{2}}+1)\left[( \sqrt a+\sqrt b)^{2F(n)}+(\sqrt a- \sqrt b)^{2F(n)}\right]$$
其中 $$F(n)= \begin{cases} 1 &\mbox{n=0,1}\\ F(n-1) + F(n-2)&\mbox{n>1} \end{cases}$$
其中 $1\leq a,b,n,p \leq 2*10^9,p是奇素数，且a,b < p$

思路：

这个式子的前两项$(a^{\frac{p-1}{2}}+1)(b^{\frac{p-1}{2}}+1)$用一个简单的快速幂就可以解决，关键在于后面的那一项
根据数列的特征方程的知识，可以将其转化为递推式，从而用矩阵快速幂来解决（关于数列的特征方程，不会的可以看一下这里）
这里稍微说一下，这里怎么得到的递推式
我们先不考虑$F(n)$，因为这是一个斐波那契数列，一会再说
设

$$g(n) = ( \sqrt a+\sqrt b)^{2n}+(\sqrt a- \sqrt b)^{2n}$$
稍微化简一下可以得到 $$g(n) = (a+b+2\sqrt {ab})^n+(a+b-\sqrt {ab})^n$$
根据 $x_1 = (a+b+2\sqrt {ab})\ \ \ x_2 = (a+b-2\sqrt {ab})$
得 $x_1+x_2 = 2(a+b) \ \ \ \ \ x_1x_2 = (a-b)^2$
构造 特征方程 $$x^2 -2(a+b)x + (a-b)^2 = 0$$
可以得到递推式
$$g(n) = 2(a+b) g_{n-1} - (a-b)^2g_{n-2}$$
构造矩阵 $$\begin{bmatrix}g(n)&g(n-1) \end{bmatrix} = \begin{bmatrix}g(n-1)&g(n-2) \end{bmatrix}\begin{bmatrix}2(a+b)&1\\-(a-b)^2&0 \end{bmatrix}$$
这里要注意一下，矩阵里面有一个负数，我们在初始化矩阵的时候就根据mod的值将其变成正数
$g(n)$这一个我们解决了，然后就是求 $g(n)的第F(n)项了，即g(F(n))$,其中 $F(n)$是斐波那契数列，数值会非常大，也要用矩阵快速幂来解决，在计算的过程中要取模，但由于F(n)是在指数上的，就不能随便取模，也就涉及到了降幂公式 $$A^{B}\equiv\ A^{B\ mod \ \varphi(C)\ +\ \varphi(C)} \ (mod\ C) \ \ \ if\ B \geq\varphi(C)$$
其中 $\varphi(C)$为欧拉函数，这里由于取模的数p是质数，不必考虑 $B \geq\varphi(C)$这个条件，而且素数p的欧拉函数值即为p-1
也就是说，我们在计算 $F(n)$的值的时候应该对 $p-1$取模，而不是p
到这里这个题就基本上解决了，然后就是实现的问题了，可以参考一下我的代码

---

具体代码如下：
Result：Accepted$\ \ \ $Memory: 1740K$\ \ \ $Time : 0MS

#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
ll A,B,n,p;
int T;
ll ans;
ll a[2][2],b[2][2];
ll mod_pow(ll x,ll n, ll mod)
{
    x%=mod;
    ll res = 1;
    while(n>0)
    {
        if(n&1) res = res*x%mod;
        x = x*x%mod;
        n>>=1;
    }
    return res;
}
void mmul(ll a[2][2],ll b[2][2],ll s[2][2],ll mod)
{
    ll tmp[2][2];
    for(int i=0;i<2;i++)
       for(int j=0;j<2;j++)
       {
          tmp[i][j]=0;
          for(int k=0;k<2;k++)
             tmp[i][j]=(tmp[i][j]+a[i][k]*b[k][j]%mod+mod)%mod;
         }
    for(int i=0;i<2;i++)
       for(int j=0;j<2;j++)
          s[i][j]=tmp[i][j];
}
ll f(ll n,ll mod)
{
    memset(a,0,sizeof(a));
    for(int i=0;i<2;i++)a[i][i]=1;
    memset(b,0,sizeof(b));
    b[0][0]=1;
    b[0][1]=1;
    b[1][0]=1;
    b[1][1]=0;
    n--;
    while(n)
    {
        if(n&1)mmul(a,b,a,mod);
        mmul(b,b,b,mod);
        n>>=1;
    }
    return (a[0][0]+a[0][1])%mod;
}
ll g(ll n,ll mod)
{
    memset(a,0,sizeof a);
    for(int i=0;i<2;i++) a[i][i]=1;
    memset(b,0,sizeof b);
    b[0][0] = 2*(A+B);
    b[0][1] = 1;
    b[1][0] = ((-(A-B)*(A-B))%mod+mod)%mod;//将矩阵中的负数处理成正数
    b[1][1] = 0;
    n--;
    n += p-1;
    while(n)
    {
        if(n&1) mmul(a,b,a,mod);
        mmul(b,b,b,mod);
        n>>=1;
    }
    return (a[0][0]*2*(A+B)%mod+2*a[1][0]%mod)%mod;
}
int main()
{
    ios::sync_with_stdio(false);
    cin>>T;
    while(T--)
    {
        cin>>A>>B>>n>>p;
        ans = (mod_pow(A,(p-1)/2,p)+1)%p;
        ans = ans*(mod_pow(B,(p-1)/2,p)+1)%p;
        ans = ans*(g(f(n,p-1),p)%p)%p;
        cout<<ans<<endl;
    }
}