HDU4565 So Easy!【矩阵连乘】【推导】

题目链接：

http://acm.hdu.edu.cn/showproblem.php?pid=4565

题目大意：

定义Sn = [( a + √b )^n] % m，[x]表示x向上取整，比如[3.14] = 4。给你a b n m的值，

求Sn是多少。

思路：

这道题很经典，因为(a-1)^2< b < a^2，所以0 < |a-√(b)| < 1，所以

Sn = [( a + √b )^n] % m = ( [( a + √b )^n]  + [(
a - √b))^n]  ) % m。

即右边其实是一个整数，如果将右边二项式展开的话，除了相互抵消的部分，剩下的部分全为

整数。这个式子设An = (a + √b)^n，Bn = (a - √b)^n，Cn = [An]，Sn = Cn % m。

先来求Cn的递推公式。

设Cn = pC(n-1) + qC(n-2)，特征方程为x^2 = p*x + q，从[( a + √b )^n]  + [( a - √b))^n]

可以看出特征根分别为a + √b和a - √b。则p = 2*a，q = b - a^2。

则Cn = 2*a*C(n-1) + (b-a^2)C(n-2)。然后开始构造矩阵

[   Cn   ]   =  [ 2*a    b-a^2 ] * [C(n-1)]

[C(n-1)]       [    1            0   ]   [C(n-2)]

[   Cn   ]   =  [ 2*a    b-a^2 ] ^(n-2)   *   [C2]

[C(n-1)]       [    1            0   ]                   [C1]

由上边Sn公式可知，C2 = 2*(a^2 + b)，C1 = 2*a。

然后用矩阵快速幂求解得到Cn，然后输出就可以了。

AC代码：

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
#include<cmath>
#define LL __int64
using namespace std;
const int MAXN = 2;
struct Matrix
{
    LL m[MAXN][MAXN];
};
LL A,B,N,M,ret,mo;
Matrix I =
{
    1,0,
    0,1
};
Matrix Multi(Matrix a, Matrix b,LL mo)
{
    Matrix c;
    for(int i = 0; i < MAXN; ++i)
    {
        for(int j = 0; j < MAXN; ++j)
        {
            c.m[i][j] = 0;
            for(int k = 0; k < MAXN; ++k)
            {
                c.m[i][j] += a.m[i][k] * b.m[k][j];
            }
            c.m[i][j] %= mo;
        }
    }
    return c;
}

Matrix Power(Matrix a, LL k)
{
    Matrix ans = I, p = a;
    while(k)
    {
        if(k&1)
            ans = Multi(ans,p,mo);
        p = Multi(p,p,mo);
        k >>= 1;
    }
    return ans;
}

int main()
{

    Matrix a,ans;
    while(cin >> A >> B >> N >> mo)
    {
        a.m[0][0] = 2*A%mo;
        a.m[0][1] = ((B%mo-A*A%mo)%mo + mo) % mo;
        a.m[1][0] = 1;
        a.m[1][1] = 0;
        if(N == 1)
            cout << 2*A % mo << endl;
        else
        {
            ans = Power(a,N-2);
            ret = (ans.m[0][0]%mo*2*(A*A%mo+B%mo)%mo + 2*A%mo*ans.m[0][1]%mo)%mo;
            ret %= mo;
            cout << ret << endl;
        }
    }

    return 0;
}