正文

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

IP云
IP云V管理员/4 阅读
11/29
提示：扫一扫查出行
【扫一扫了解最新限行尾号】
复制提示

题目: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 $使用矩阵快速幂可以快速计算。 

 #include <iostream>

 #include <cstdio>

 using namespace std;

 const int maxN=;

 int MOD;

 struct Matrix{int n,a[maxN][maxN];}A,E;

 Matrix operator+(Matrix A,Matrix B){

     Matrix RT{};RT.n=A.n;

     for(int i=;i<RT.n;++i)

     for(int j=;j<RT.n;++j)

         RT.a[i][j]=(A.a[i][j]+B.a[i][j])%MOD;

     return RT;

 }

 Matrix operator*(Matrix A,Matrix B){

     Matrix RT{};RT.n=A.n;

     for(int i=;i<A.n;++i)

     for(int j=;j<A.n;++j)

         for(int k=;k<A.n;++k)

             (RT.a[i][j]+=A.a[i][k]*B.a[k][j])%=MOD;

     return RT;

 }

 Matrix operator^(Matrix A,int n){

     Matrix RT=E;

     for(;n;n>>=){

         if(n&)RT=RT*A;

         A=A*A;

     }

     return RT;

 }

 Matrix Sum(Matrix &A,int n){

     if(n==)return A;

     if(n&)return (A^n) + Sum(A,n-);

     Matrix B=Sum(A,n>>);

     return B+((A^(n>>))*B);

 }

 int main(){

     for(int n,k;~scanf("%d %d %d",&n,&k,&MOD);){

         A.n=E.n=n;

         for(int i=;i<n;++i)

         for(int j=;j<n;++j){

             scanf("%d",&A.a[i][j]);

             A.a[i][j]%=MOD;E.a[i][j]=;

         }

         for(int i=;i<n;++i)E.a[i][i]=;

         Matrix RT=Sum(A,k);

         for(int i=;i<n;++i){

             for(int j=;j<n;++j)

                 printf("%d ",RT.a[i][j]);

             puts("");

         }

     }

     return ;

 }

方法二:对于多项式，有秦九韶算法，而对$ 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] $ 

 #include <iostream>

 #include <cstdio>

 using namespace std;

 const int maxN=;

 int MOD;

 struct Matrix{int n,a[maxN][maxN];}A[],B[];

 Matrix operator+(Matrix A,Matrix B){

     Matrix RT{};RT.n=A.n;

     for(int i=;i<RT.n;++i)

     for(int j=;j<RT.n;++j)

         RT.a[i][j]=(A.a[i][j]+B.a[i][j])%MOD;

     return RT;

 }

 Matrix operator*(Matrix A,Matrix B){

     Matrix RT{};RT.n=A.n;

     for(int i=;i<A.n;++i)

     for(int j=;j<A.n;++j)

         for(int k=;k<A.n;++k)

             (RT.a[i][j]+=A.a[i][k]*B.a[k][j])%=MOD;

     return RT;

 }

 int main(){

     Matrix E{};

     for(int n,k;~scanf("%d %d %d",&n,&k,&MOD);){

         A[].n=E.n=n;

         for(int i=;i<n;++i)

         for(int j=;j<n;++j){

             scanf("%d",&A[].a[i][j]);

             A[].a[i][j]%=MOD;E.a[i][j]=;

         }

         for(int i=;i<n;++i)E.a[i][i]=;

         for(int i=;k>>i;++i)A[i]=A[i-]*A[i-];

         B[]=A[];

         for(int i=;k>>i;++i)B[i]=B[i-]*(A[i-]+E);

         Matrix RT{};RT.n=n;

         for(int i=;~i;--i)if(k>>i&)RT=RT*A[i]+B[i];

         for(int i=;i<n;++i){

             for(int j=;j<n;++j)

                 printf("%d ",RT.a[i][j]);

             puts("");

         }

     }

     return ;

 }

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

 #include <iostream>

 #include <cstdio>

 using namespace std;

 const int maxN=;

 int MOD;

 struct Matrix{int n,a[maxN][maxN];}A[],B[];

 Matrix operator+(Matrix A,Matrix B){

     Matrix RT{};RT.n=A.n;

     for(int i=;i<RT.n;++i)

     for(int j=;j<RT.n;++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{};RT.n=A.n;

     for(int i=;i<A.n;++i)

     for(int j=;j<A.n;++j){

         for(int k=;k<A.n;++k)

             RT.a[i][j]+=A.a[i][k]*B.a[k][j];

         RT.a[i][j]%=MOD;

     }

     return RT;

 }

 int main(){

     Matrix E{};

     for(int n,k;~scanf("%d %d %d",&n,&k,&MOD);){

         A[].n=E.n=n;

         for(int i=;i<n;++i)

         for(int j=;j<n;++j){

             scanf("%d",&A[].a[i][j]);

             A[].a[i][j]%=MOD;E.a[i][j]=;

         }

         for(int i=;i<n;++i)E.a[i][i]=;

         for(int i=;k>>i;++i)A[i]=A[i-]*A[i-];

         B[]=A[];

         for(int i=;k>>i;++i)B[i]=B[i-]*(A[i-]+E);

         Matrix RT{};RT.n=n;

         for(int i=;~i;--i)if(k>>i&)RT=RT*A[i]+B[i];

         for(int i=;i<n;++i){

             for(int j=;j<n;++j)

                 printf("%d ",RT.a[i][j]);

             puts("");

         }

     }

     return ;

 }

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

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

 #include <iostream>

 #include <cstdio>

 using namespace std;

 const int maxN=;

 int MOD;

 struct Matrix{int n,a[maxN][maxN];}A;

 Matrix operator+(Matrix A,Matrix B){

     Matrix RT{};RT.n=A.n;

     for(int i=;i<RT.n;++i)

     for(int j=;j<RT.n;++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{};RT.n=A.n;

     for(int i=;i<A.n;++i)

     for(int j=;j<A.n;++j){

         for(int k=;k<A.n;++k)

             RT.a[i][j]+=A.a[i][k]*B.a[k][j];

         RT.a[i][j]%=MOD;

     }

     return RT;

 }

 Matrix Sum(Matrix A,int k){

     Matrix RT{};int n=RT.n=A.n;

     for(int t=n/,i=;i<n;++i)RT.a[i][i+t]=;

     for(;k;k>>=){

         if(k&)RT=RT*A;

         A=A*A;

     }

     return RT;

 }

 int main(){

     for(int n,k;~scanf("%d %d %d",&n,&k,&MOD);){

         A.n=n<<;

         for(int i=;i<n;++i){

             A.a[i][i]=;

             for(int j=,t;j<n;++j){

                 scanf("%d",&t);if(t>=MOD)t%=MOD;

                 A.a[i+n][j]=A.a[i+n][j+n]=t;

             }

         }

         Matrix RT=Sum(A,k);

         for(int i=;i<n;++i){

             for(int j=;j<n;++j)

                 printf("%d ",RT.a[i][j]);

             puts("");

         }

     }

     return ;

 }

题目: 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) $

就和上面一样了。 

 #include <iostream>

 #include <cstdio>

 using namespace std;

 typedef long long LL;

 const int maxN=;

 LL MOD;

 struct Matrix{LL a[maxN][maxN];}A[],B[],E,Fib;

 Matrix operator+(Matrix A,Matrix B){

     Matrix RT{};

     for(int i=;i<;++i)

     for(int j=;j<;++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{};

     for(int i=;i<;++i)

     for(int j=;j<;++j){

         for(int k=;k<;++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=E;

     for(;n;n>>=){

         if(n&)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.a[][]=;Fib.a[][]=;

         Fib.a[][]=;Fib.a[][]=;

         E.a[][]=E.a[][]=;

         E.a[][]=E.a[][]=;

         A[]=Fib^k;

         for(int i=;n>>i;++i)A[i]=A[i-]*A[i-];

         B[]=A[];

         for(int i=;n>>i;++i)B[i]=B[i-]*(A[i-]+E);

         Matrix RT{};

         for(int i=;~i;--i)if(n>>i&)RT=RT*A[i]+B[i];

         RT=(RT+E)*(Fib^b);

         printf("%lld\n",RT.a[][]);

     }

     return ;

 }
#include <iostream>

#include <cstdio>

using namespace std;

typedef long long LL;

const int maxN=;

LL MOD;

struct Matrix{int n;LL a[maxN][maxN];}A,Fib;

Matrix operator+(Matrix A,Matrix B){

    Matrix RT{};

    for(int i=;i<;++i)

    for(int j=;j<;++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{};int n=RT.n=A.n;

    for(int i=;i<n;++i)

    for(int j=;j<n;++j){

        for(int k=;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=;i<RT.n;++i)for(int j=;j<RT.n;++j)RT.a[i][j]=;

    for(int i=;i<RT.n;++i)RT.a[i][i]=;

    for(;n;n>>=){

        if(n&)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=;

        Fib.a[][]=;Fib.a[][]=;

        Fib.a[][]=;Fib.a[][]=;

        A=Fib^k;

        A.n=;

        A.a[][]=A.a[][]=A.a[][];

        A.a[][]=A.a[][]=A.a[][];

        A.a[][]=A.a[][]=A.a[][];

        A.a[][]=A.a[][]=A.a[][];

        A.a[][]=;A.a[][]=;

        A.a[][]=;A.a[][]=;

        A=A^n;

        A.n=;

        A.a[][]=A.a[][]+;A.a[][]=A.a[][]+;

        A.a[][]=A.a[][]+;A.a[][]=A.a[][]+;

        A=(Fib^b)*A;

        printf("%lld\n",A.a[][]);

    }

    return ;

}

POJ3233]Matrix Power Series && [HDU1588]Gauss Fibonacci的更多相关文章
  1. [POJ3233]Matrix Power Series 分治+矩阵

    本文为博主原创文章,欢迎转载,请注明出处 www.cnblogs.com/yangyaojia [POJ3233]Matrix Power Series 分治+矩阵 题目大意 A为n×n(n<= ...

  2. POJ3233 Matrix Power Series 矩阵快速幂 矩阵中的矩阵

    Matrix Power Series Time Limit: 3000MS   Memory Limit: 131072K Total Submissions: 27277   Accepted:  ...

  3. POJ3233 Matrix Power Series

    Description Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak. ...

  4. POJ3233:Matrix Power Series(矩阵快速幂+二分）

    http://poj.org/problem?id=3233 题目大意:给定矩阵A,求A + A^2 + A^3 + … + A^k的结果(两个矩阵相加就是对应位置分别相加).输出的数据mod m.k ...

  5. POJ3233 Matrix Power Series(矩阵快速幂+分治)

    Description Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak. ...

  6. POJ3233:Matrix Power Series(矩阵快速幂+递推式)

    传送门 题意 给出n,m,k,求 \[\sum_{i=1}^kA^i\] A是矩阵 分析 我们首先会想到等比公式,然后得到这样一个式子: \[\frac{A^{k+1}-E}{A-E}\] 发现要用矩 ...

  7. POJ-3233 Matrix Power Series 矩阵A^1+A^2+A^3...求和转化

    S(k)=A^1+A^2...+A^k. 保利求解就超时了,我们考虑一下当k为偶数的情况,A^1+A^2+A^3+A^4...+A^k,取其中前一半A^1+A^2...A^k/2,后一半提取公共矩阵A ...

  8. POJ3233 Matrix Power Series（快速幂求等比矩阵和）

    题面 \(solution:\) 首先,如果题目只要我们求\(A^K\) 那这一题我们可以直接模版矩乘快速幂来做,但是它现在让我们求$\sum_{i=1}^{k}{(A^i)} $ 所以我们思考一下这 ...

  9. poj3233 Matrix Power Series（矩阵快速幂）

    题目要求的是 A+A2+...+Ak,而不是单个矩阵的幂. 那么可以构造一个分块的辅助矩阵 S,其中 A 为原矩阵,E 为单位矩阵,O 为0矩阵    将 S 取幂,会发现一个特性: Sk +1右上角 ...

随机推荐
  1. mariadb（三）查

    -查询基本使用(条件,排序,聚合函数,分组,分页) 1)创建一个表结构然后添加数据 create table baba (id int unsigned not null auto_increment ...

  2. curl 和 wget 命令

    1. curl curl 支持 HTTP.HTTPS.FTP 等协议,还支持 POST.cookies.认证.从指定偏移处下载部分文件.User-Agent.限速.文件大小.进度条等特征. 1.1 选 ...

  3. 使用定时器settimeout、setInterval执行能传递参数的函数

    无论是window.setTimeout还是window.setInterval,在使用函数名作为调用句柄时都不能带参数,而在许多场合必须要带参数,这就需要想方法解决.经网上查询后整理如下: 例如对于 ...

  4. 舔狗【2019河北省大学生程序设计竞赛 J题】

    题目描述 > “舔狗舔狗,> 舔到最后,> 一无所有.” 有 n 只舔狗,每只舔狗的心中都有自己朝思暮想的一位. 每个人虽然受到了一万次拒绝,还毅然第一万零一次鼓起勇气. 作为一个不 ...

  5. Django+python报错与解决方法

    编码信息导致修改django文件下的debug.py文件331行的编码格式 with Path(CURRENT_DIR, 'templates', 'technical_500.html').open ...

  6. 《JAVA设计模式》之模板模式(Template)

    在阎宏博士的<JAVA与模式>一书中开头是这样描述模板方法(Template Method)模式的: 模板方法模式是类的行为模式.准备一个抽象类,将部分逻辑以具体方法以及具体构造函数的形式 ...

  7. vscode内置c++ define选项

    https://stackoverflow.com/questions/46415076/how-can-i-define-macros-for-the-c-intellisense-engine

  8. 2.openshift授权策略和scc理解

    policy管理概念 查看policy的方式可以通过cli进行查看 Roles grant various levels of access in the system-wide cluster po ...

  9. bzoj1779 [Usaco2010 Hol]Cowwar 奶牛战争（网络流）

    1779: [Usaco2010 Hol]Cowwar 奶牛战争 Time Limit: 10 Sec  Memory Limit: 64 MBSubmit: 302  Solved: 131[Sub ...

  10. C# 跨线程调用控件的4中方法

    原文:C# 跨线程调用控件 在C# 的应用程序开发中, 我们经常要把UI线程和工作线程分开,防止界面停止响应.  同时我们又需要在工作线程中更新UI界面上的控件, 下面介绍几种常用的方法 阅读目录 线 ...

相关阅读

目录[+]