POJ3233]Matrix Power Series && [HDU1588]Gauss Fibonacci

题目:Matrix Power Series

传送门：http://poj.org/problem?id=3233

分析：

方法一：引用Matrix67大佬的矩阵十题：这道题两次二分，相当经典。首先我们知道，A^i可以二分求出。然后我们需要对整个题目的数据规模k进行二分。比如，当k=6时，有：$ S(6)= A + A^2 + A^3 + A^4 + A^5 + A^6 =\underline{(A + A^2 + A^3)} + A^3*\underline{(A + A^2 + A^3)}。   $

即对于k：如果k是偶数，就二分减小规模，$ S(k)=S(\frac{k}{2})+A^{\frac{n}{2}} *S(\frac{k}{2}) $。如果k是奇数，那么 $ S(k)=S(k-1)+A^n $; 其中 $ A^n $使用矩阵快速幂可以快速计算。

 1 #include <iostream>
 2 #include <cstdio>
 3 using namespace std;
 4 const int maxN=35;
 5 int MOD;
 6 struct Matrix{int n,a[maxN][maxN];}A,E;
 7 Matrix operator+(Matrix A,Matrix B){
 8     Matrix RT{0};RT.n=A.n;
 9     for(int i=0;i<RT.n;++i)
10     for(int j=0;j<RT.n;++j)
11         RT.a[i][j]=(A.a[i][j]+B.a[i][j])%MOD;
12     return RT;
13 }
14 Matrix operator*(Matrix A,Matrix B){
15     Matrix RT{0};RT.n=A.n;
16     for(int i=0;i<A.n;++i)
17     for(int j=0;j<A.n;++j)
18         for(int k=0;k<A.n;++k)
19             (RT.a[i][j]+=A.a[i][k]*B.a[k][j])%=MOD;
20     return RT;
21 }
22 Matrix operator^(Matrix A,int n){
23     Matrix RT=E;
24     for(;n;n>>=1){
25         if(n&1)RT=RT*A;
26         A=A*A;
27     }
28     return RT;
29 }
30 Matrix Sum(Matrix &A,int n){
31     if(n==1)return A;
32     if(n&1)return (A^n) + Sum(A,n-1);
33     Matrix B=Sum(A,n>>1);
34     return B+((A^(n>>1))*B);
35 }
36 int main(){
37     for(int n,k;~scanf("%d %d %d",&n,&k,&MOD);){
38         A.n=E.n=n;
39         for(int i=0;i<n;++i)
40         for(int j=0;j<n;++j){
41             scanf("%d",&A.a[i][j]);
42             A.a[i][j]%=MOD;E.a[i][j]=0;
43         }
44         for(int i=0;i<n;++i)E.a[i][i]=1;
45         Matrix RT=Sum(A,k);
46         for(int i=0;i<n;++i){
47             for(int j=0;j<n;++j)
48                 printf("%d ",RT.a[i][j]);
49             puts("");
50         }
51     }
52     return 0;
53 }
View Code

方法二:对于多项式，有秦九韶算法，而对$ A + A^2 + A^3 + … + A^k $ 这样的多项式，可以二进制优化秦九韶算法。$ S( 2^i ) $ 是可以快速计算的: $S( 2^i ) = S( 2^{(i-1)} ) * (E + A^{(2^i)} )$。记：$A[i]=A^{2^i}$、$S[i]=S(2^i)$;预处理A[i]、s[i]。

比如k=6时：$ S(6)=\underline{(A + A^2 + A^3 + A^ 4)} * A^2 + \underline{(A + A^2 )} = S(4) * A(2) + S(2) = S[2]*A[1]+S[1] $

 1 #include <iostream>
 2 #include <cstdio>
 3 using namespace std;
 4 const int maxN=35;
 5 int MOD;
 6 struct Matrix{int n,a[maxN][maxN];}A[35],B[35];
 7 Matrix operator+(Matrix A,Matrix B){
 8     Matrix RT{0};RT.n=A.n;
 9     for(int i=0;i<RT.n;++i)
10     for(int j=0;j<RT.n;++j)
11         RT.a[i][j]=(A.a[i][j]+B.a[i][j])%MOD;
12     return RT;
13 }
14 Matrix operator*(Matrix A,Matrix B){
15     Matrix RT{0};RT.n=A.n;
16     for(int i=0;i<A.n;++i)
17     for(int j=0;j<A.n;++j)
18         for(int k=0;k<A.n;++k)
19             (RT.a[i][j]+=A.a[i][k]*B.a[k][j])%=MOD;
20     return RT;
21 }
22 int main(){
23     Matrix E{0};
24     for(int n,k;~scanf("%d %d %d",&n,&k,&MOD);){
25         A[0].n=E.n=n;
26         for(int i=0;i<n;++i)
27         for(int j=0;j<n;++j){
28             scanf("%d",&A[0].a[i][j]);
29             A[0].a[i][j]%=MOD;E.a[i][j]=0;
30         }
31         for(int i=0;i<n;++i)E.a[i][i]=1;
32         for(int i=1;k>>i;++i)A[i]=A[i-1]*A[i-1];
33         B[0]=A[0];
34         for(int i=1;k>>i;++i)B[i]=B[i-1]*(A[i-1]+E);
35         Matrix RT{0};RT.n=n;
36         for(int i=30;~i;--i)if(k>>i&1)RT=RT*A[i]+B[i];
37         for(int i=0;i<n;++i){
38             for(int j=0;j<n;++j)
39                 printf("%d ",RT.a[i][j]);
40             puts("");
41         }
42     }
43     return 0;
44 }
View Code

优化一下常数（141s -> 32s)：

 1 #include <iostream>
 2 #include <cstdio>
 3 using namespace std;
 4 const int maxN=35;
 5 int MOD;
 6 struct Matrix{int n,a[maxN][maxN];}A[35],B[35];
 7 Matrix operator+(Matrix A,Matrix B){
 8     Matrix RT{0};RT.n=A.n;
 9     for(int i=0;i<RT.n;++i)
10     for(int j=0;j<RT.n;++j){
11         RT.a[i][j]=A.a[i][j]+B.a[i][j];
12         if(RT.a[i][j]>=MOD)RT.a[i][j]-=MOD;
13     }
14     return RT;
15 }
16 Matrix operator*(Matrix A,Matrix B){
17     Matrix RT{0};RT.n=A.n;
18     for(int i=0;i<A.n;++i)
19     for(int j=0;j<A.n;++j){
20         for(int k=0;k<A.n;++k)
21             RT.a[i][j]+=A.a[i][k]*B.a[k][j];
22         RT.a[i][j]%=MOD;
23     }
24     return RT;
25 }
26 int main(){
27     Matrix E{0};
28     for(int n,k;~scanf("%d %d %d",&n,&k,&MOD);){
29         A[0].n=E.n=n;
30         for(int i=0;i<n;++i)
31         for(int j=0;j<n;++j){
32             scanf("%d",&A[0].a[i][j]);
33             A[0].a[i][j]%=MOD;E.a[i][j]=0;
34         }
35         for(int i=0;i<n;++i)E.a[i][i]=1;
36         for(int i=1;k>>i;++i)A[i]=A[i-1]*A[i-1];
37         B[0]=A[0];
38         for(int i=1;k>>i;++i)B[i]=B[i-1]*(A[i-1]+E);
39         Matrix RT{0};RT.n=n;
40         for(int i=30;~i;--i)if(k>>i&1)RT=RT*A[i]+B[i];
41         for(int i=0;i<n;++i){
42             for(int j=0;j<n;++j)
43                 printf("%d ",RT.a[i][j]);
44             puts("");
45         }
46     }
47     return 0;
48 }

 

方法三：有种很有意思的打法：构造矩阵

$ T = \left[ \begin{array}{cc} E & 0 \\ A & A \end{array}  \right]$ 、 $ T^2 =  \left[ \begin{array}{cc} E & 0 \\ {A + A^2}  & A^2 \end{array}  \right]$ 、 $ T^2 =  \left[ \begin{array}{cc} E & 0 \\ {A + A^2 + A^3}  & A^3 \end{array} \right]$ 、

 这个矩阵的左下角就是矩阵次方和。$ T^n =  \left[ \begin{array}{cc} E & 0 \\ {A + A^2 + .. A^n}  & A^n \end{array} \right]$

这个矩阵可以直接快速幂求。

把RT矩阵设为 $ RT =  \left[ \begin{array}{cc} 0 & E \\ 0  & 0 \end{array} \right] $，然后 $ RT \times T^n =   \left[ \begin{array}{cc} {A + A^2 + .. A^n} & 0 \\ 0  & 0 \end{array} \right] $

 1 #include <iostream>
 2 #include <cstdio>
 3 using namespace std;
 4 const int maxN=65;
 5 int MOD;
 6 struct Matrix{int n,a[maxN][maxN];}A;
 7 Matrix operator+(Matrix A,Matrix B){
 8     Matrix RT{0};RT.n=A.n;
 9     for(int i=0;i<RT.n;++i)
10     for(int j=0;j<RT.n;++j){
11         RT.a[i][j]=A.a[i][j]+B.a[i][j];
12         if(RT.a[i][j]>=MOD)RT.a[i][j]-=MOD;
13     }
14     return RT;
15 }
16 Matrix operator*(Matrix A,Matrix B){
17     Matrix RT{0};RT.n=A.n;
18     for(int i=0;i<A.n;++i)
19     for(int j=0;j<A.n;++j){
20         for(int k=0;k<A.n;++k)
21             RT.a[i][j]+=A.a[i][k]*B.a[k][j];
22         RT.a[i][j]%=MOD;
23     }
24     return RT;
25 }
26 Matrix Sum(Matrix A,int k){
27     Matrix RT{0};int n=RT.n=A.n;
28     for(int t=n/2,i=0;i<n;++i)RT.a[i][i+t]=1;
29     for(;k;k>>=1){
30         if(k&1)RT=RT*A;
31         A=A*A;
32     }
33     return RT;
34 }
35 int main(){
36     for(int n,k;~scanf("%d %d %d",&n,&k,&MOD);){
37         A.n=n<<1;
38         for(int i=0;i<n;++i){
39             A.a[i][i]=1;
40             for(int j=0,t;j<n;++j){
41                 scanf("%d",&t);if(t>=MOD)t%=MOD;
42                 A.a[i+n][j]=A.a[i+n][j+n]=t;
43             }
44         }
45         Matrix RT=Sum(A,k);
46         for(int i=0;i<n;++i){
47             for(int j=0;j<n;++j)
48                 printf("%d ",RT.a[i][j]);
49             puts("");
50         }
51     }
52     return 0;
53 }

题目: Gauss Fibonacci

链接:http://acm.hdu.edu.cn/showproblem.php?pid=1588

分析:Fibonacci数列可用矩阵:$ Fib =  \left[ \begin{array}{cc} 0 & 1 \\ 1  & 1 \end{array} \right] $ 表示,f(i)=$Fib^n$的左下角元素的值。

然后即求 $ Fib^b + Fib^{k+b} + Fib^{2*k+b} + ... Fib^{(n-1)*k+b} $ 

即:$ Fib^b * (E + Fib^k + Fib^{2*k} + Fib^{3*k} + ... +Fib^{(n-1)*k} $ = $ Fib^b * (E + (Fib^k) + {(Fib^k)}^2 + {(Fib^k)}^3 + ... +{(Fib^k)}^{(n-1)} )$

令  $ A=Fib^k $;

则求: $ Fib^b * (E + A + A^2 + A^3 +.. A^n) $

就和上面一样了。

 

 1 #include <iostream>
 2 #include <cstdio>
 3 using namespace std;
 4 typedef long long LL;
 5 const int maxN=2;
 6 LL MOD;
 7 struct Matrix{LL a[maxN][maxN];}A[35],B[35],E,Fib;
 8 Matrix operator+(Matrix A,Matrix B){
 9     Matrix RT{0};
10     for(int i=0;i<2;++i)
11     for(int j=0;j<2;++j){
12         RT.a[i][j]=A.a[i][j]+B.a[i][j];
13         if(RT.a[i][j]>=MOD)RT.a[i][j]-=MOD;
14     }
15     return RT;
16 }
17 Matrix operator*(Matrix A,Matrix B){
18     Matrix RT{0};
19     for(int i=0;i<2;++i)
20     for(int j=0;j<2;++j){
21         for(int k=0;k<2;++k)
22             RT.a[i][j]+=A.a[i][k]*B.a[k][j];
23         RT.a[i][j]%=MOD;
24     }
25     return RT;
26 }
27 Matrix operator^(Matrix A,LL n){
28     Matrix RT=E;
29     for(;n;n>>=1){
30         if(n&1)RT=RT*A;
31         A=A*A;
32     }
33     return RT;
34 }
35 int main(){
36     for(LL k,b,n;~scanf("%lld %lld %lld %lld",&k,&b,&n,&MOD);){
37         --n;
38         Fib.a[0][0]=0;Fib.a[0][1]=1;
39         Fib.a[1][0]=1;Fib.a[1][1]=1;
40         E.a[0][0]=E.a[1][1]=1;
41         E.a[1][0]=E.a[0][1]=0;
42 
43         A[0]=Fib^k;
44         for(int i=1;n>>i;++i)A[i]=A[i-1]*A[i-1];
45         B[0]=A[0];
46         for(int i=1;n>>i;++i)B[i]=B[i-1]*(A[i-1]+E);
47         
48         Matrix RT{0};
49         for(int i=30;~i;--i)if(n>>i&1)RT=RT*A[i]+B[i];
50         RT=(RT+E)*(Fib^b);
51         printf("%lld\n",RT.a[1][0]);
52     }
53     return 0;
54 }
View Code

 

#include <iostream>
#include <cstdio>
using namespace std;
typedef long long LL;
const int maxN=4;
LL MOD;
struct Matrix{int n;LL a[maxN][maxN];}A,Fib;
Matrix operator+(Matrix A,Matrix B){
    Matrix RT{0};
    for(int i=0;i<2;++i)
    for(int j=0;j<2;++j){
        RT.a[i][j]=A.a[i][j]+B.a[i][j];
        if(RT.a[i][j]>=MOD)RT.a[i][j]-=MOD;
    }
    return RT;
}
Matrix operator*(Matrix A,Matrix B){
    Matrix RT{0};int n=RT.n=A.n;
    for(int i=0;i<n;++i)
    for(int j=0;j<n;++j){
        for(int k=0;k<n;++k)
            RT.a[i][j]+=A.a[i][k]*B.a[k][j];
        RT.a[i][j]%=MOD;
    }
    return RT;
}
Matrix operator^(Matrix A,LL n){
    Matrix RT;
    RT.n=A.n;
    for(int i=0;i<RT.n;++i)for(int j=0;j<RT.n;++j)RT.a[i][j]=0;
    for(int i=0;i<RT.n;++i)RT.a[i][i]=1;
    for(;n;n>>=1){
        if(n&1)RT=RT*A;
        A=A*A;
    }
    return RT;
}
int main(){
    for(LL k,b,n;~scanf("%lld %lld %lld %lld",&k,&b,&n,&MOD);){
        --n;
        Fib.n=2;
        Fib.a[0][0]=0;Fib.a[0][1]=1;
        Fib.a[1][0]=1;Fib.a[1][1]=1;
        
        A=Fib^k;
        A.n=4;
        A.a[2][0]=A.a[2][2]=A.a[0][0];
        A.a[2][1]=A.a[2][3]=A.a[0][1];
        A.a[3][0]=A.a[3][2]=A.a[1][0];
        A.a[3][1]=A.a[3][3]=A.a[1][1];
        A.a[0][0]=1;A.a[0][1]=0;
        A.a[1][0]=0;A.a[1][1]=1;

        A=A^n;
        A.n=2;
        A.a[0][0]=A.a[2][0]+1;A.a[0][1]=A.a[2][1]+0;
        A.a[1][0]=A.a[3][0]+0;A.a[1][1]=A.a[3][1]+1;

        A=(Fib^b)*A;
        printf("%lld\n",A.a[1][0]);
    }
    return 0;
}

 

posted @ 2018-11-18 19:48  hjj1871984569  阅读(201)  评论(0)  编辑  收藏  举报