HDU4565 So Easy! [matrix multiplication] [derivation]

Topic link ：

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

The main idea of the topic ：

Definition Sn = [( a + √b )^n] % m,[x] Express x Rounding up , such as [3.14] = 4. Here you are. a b n m Value ,

seek Sn How much is the .

Ideas ：

This is a classic question , because (a-1)^2< b < a^2, therefore 0 < |a-√(b)| < 1, therefore ^​

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

That is, the right is actually an integer , If you expand the binomial on the right , Except for the offsetting part , The rest is all

Integers . This formula sets An = (a + √b)^n,Bn = (a - √b)^n,Cn = [An],Sn = Cn % m.

First ask Cn The recurrence formula of .

set up Cn = pC(n-1) + qC(n-2), The characteristic equation is x^2 = p*x + q, from [( a + √b )^n]  + [( a - √b))^n] 

It can be seen that the characteristic roots are a + √b and a - √b. be p = 2*a,q = b - a^2.

be Cn = 2*a*C(n-1) + (b-a^2)C(n-2). Then start to construct the matrix

[   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]

From top Sn The formula shows ,C2 = 2*(a^2 + b),C1 = 2*a.

Then we use the fast power of matrix to get Cn, Then the output is OK .

AC Code ：

Log in and copy

      
      

       
       #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.