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;
       
       
}
      
      
      
      

       
       
1.
       
       
2.
       
       
3.
       
       
4.
       
       
5.
       
       
6.
       
       
7.
       
       
8.
       
       
9.
       
       
10.
       
       
11.
       
       
12.
       
       
13.
       
       
14.
       
       
15.
       
       
16.
       
       
17.
       
       
18.
       
       
19.
       
       
20.
       
       
21.
       
       
22.
       
       
23.
       
       
24.
       
       
25.
       
       
26.
       
       
27.
       
       
28.
       
       
29.
       
       
30.
       
       
31.
       
       
32.
       
       
33.
       
       
34.
       
       
35.
       
       
36.
       
       
37.
       
       
38.
       
       
39.
       
       
40.
       
       
41.
       
       
42.
       
       
43.
       
       
44.
       
       
45.
       
       
46.
       
       
47.
       
       
48.
       
       
49.
       
       
50.
       
       
51.
       
       
52.
       
       
53.
       
       
54.
       
       
55.
       
       
56.
       
       
57.
       
       
58.
       
       
59.
       
       
60.
       
       
61.
       
       
62.
       
       
63.
       
       
64.
       
       
65.
       
       
66.
       
       
67.
       
       
68.
       
       
69.
       
       
70.
       
       
71.
       
       
72.