HDU 4565 So Easy!(共轭复数 + 矩阵快速幂 数论)

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=4565

Problem Description

　　A sequence S _n is defined as:

Where a, b, n, m are positive integers.┌x┐is the ceil of x. For example, ┌3.14┐=4. You are to calculate S _n.
　　You, a top coder, say: So easy! 

 

Input

　　There are several test cases, each test case in one line contains four positive integers: a, b, n, m. Where 0< a, m < 2 ¹⁵, (a-1) ²< b < a ², 0 < b, n < 2 ³¹.The input will finish with the end of file.

 

Output

　　For each the case, output an integer S _n.

 

Sample Input

  
  

   
   2 3 1 2013
2 3 2 2013
2 2 1 2013
  
  

 

Sample Output

  
  

   
   4
14
4
  
  

 

Source

2013 ACM-ICPC长沙赛区全国邀请赛——题目重现

PS:

这类题都是通过矩阵快速幂来解决！

利用共轭复数！

学习：http://m.blog.csdn.net/blog/u011481752/26291179

讲得很好！

题意很清晰，公式直接给出了。。。经典题目了，构造共轭矩阵来做。。。

记An=(a+sqrt(b))^n,Bn=(a-sqrt(b))^n Cn=An+Bn,很容易知道。。 Cn是整数，

然后由于a b的范围可以知道Bn在0-1开区间之间，因此Sn=(Cn)%mod，接下来对Cn进行配项。。。

Cn+1=(a+sqrt(b))^n*(a+sqrt(b))+(a-sqrt(b))^n*(a-sqrt(b)),貌似没啥用，但是要充分利用共轭这个式子。

要想办法让（a+sqrt(b)）*(a-sqrt(b))=a*a-b派上用场，然后继续去把它拿去配项。。。

一开始我拿Cn*(a*a-b)未果，于是用Cn-1*(a*a-b)结合上述的Cn+1就得到了：

Cn+1=(a+sqrt(b))^n*(a+sqrt(b))+(a-sqrt(b))^n*(a-sqrt(b))

Cn-1*(a*a-b)=(a+sqrt(b))^n*(a-sqrt(b))+(a-sqrt(b))^n*(a+sqrt(b))

加减消去根号项得到：Cn+1=2*a*Cn+(b-a*a)*Cn-1，公式化简完毕，2*2阶矩阵可以直接写了：

|2*a    b-a*a|   | Cn  |  | Cn+1|
|1        0  |*  | Cn-1|= | Cn  |

因为模之后存在负数的情况，所以在MOD 的时候加上个MOD防止负数的出现。

还有两题HDU5451 and  HDU 2256 和这题类似！今年的沈阳网络赛也出现了！但是那题需要先找循环节！

代码如下：

/*
An=(a+sqrt(b))^n,Bn=(a-sqrt(b))^n   Cn=An+Bn
Cn+1=(a+sqrt(b))^n*(a+sqrt(b))+(a-sqrt(b))^n*(a-sqrt(b))
想办法让（a+sqrt(b)）*(a-sqrt(b))=a*a-b派上用场
Cn-1*(a*a-b)=(a+sqrt(b))^n*(a-sqrt(b))+(a-sqrt(b))^n*(a+sqrt(b))
*/
/*
C(n+1) = 2*a*Cn + (b-a*a)*C(n-1)，

|2*a    b-a*a|   | Cn  |  | Cn+1|
|1        0  |*  | Cn-1|= | Cn  |
*/

#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL __int64
struct Matrix
{
    LL  m[5][5];
} I, A, T;

LL a, b, n, mod;
int ssize = 2;

Matrix Mul(Matrix a,Matrix b)
{
    int i, j, k;
    Matrix c;
    for(i = 1; i <= ssize; i++)
    {
        for(j = 1; j <= ssize; j++)
        {
            c.m[i][j]=0;
            for(k = 1; k <= ssize; k++)
            {
                c.m[i][j]+=(a.m[i][k]*b.m[k][j]);
                c.m[i][j]%=mod;
            }
        }
    }
    return c;
}

Matrix quickpagow(LL n)
{
    Matrix m = A, b = I;
    while(n)
    {
        if(n & 1)
            b = Mul(b,m);
        n = n >> 1;
        m = Mul(m,m);
    }
    return b;
}
int main()
{
    LL c;
    LL s1, s2;
    while(~scanf("%I64d%I64d%I64d%I64d",&a,&b,&n,&mod))
    {
        memset(I.m,0,sizeof(I.m));
        memset(A.m,0,sizeof(A.m));
        for(int i = 1; i <= ssize; i++)
        {
            //单位矩阵
            I.m[i][i] = 1;
        }
        s1 = (2*a)%mod;
        double tt = (a+sqrt((double)b)) * (a+sqrt((double)b));
        s2 = (LL)(ceil)(tt)%mod;
        A.m[1][1] = 2*a;//初始化等比矩阵
        A.m[1][2] = b-a*a;
        A.m[2][1] = 1;

        if(n == 1)
        {
            printf("%I64d\n",s1);
            continue;
        }
        if(n == 2)
        {
            printf("%I64d\n",s2);
            continue;
        }
        T = quickpagow(n-2); //注意n-2为负的情况

        LL ans = ((T.m[1][1]*s2%mod)+mod)%mod + ((T.m[1][2]*s1%mod)+mod)%mod;

        printf("%I64d\n",ans%mod);
    }
    return 0;
}