POJ3233 Matrix Power Series 矩阵快速幂+二分求和

在看矩阵快速幂求和之前，我们先来看一下等比数列Sn=(a+a^2+a^3+...+a^n)mod M的求和取模：

实现代码如下：

#include <iostream>
#include <string.h>
#include <stdio.h>
using namespace std;
typedef long long LL;
const int M=10000009;
LL power(LL a,LL b)
{
    LL ans = 1;
    a %= M;
    while(b)
    {
        if(b & 1)
        {
            ans = ans * a % M;
            b--;
        }
        b >>= 1;
        a = a * a % M;
    }
    return ans;
}

LL sum(LL a,LL n)
{
    if(n == 1) return a;
    LL t = sum(a,n/2);
    if(n & 1)
    {
        LL cur = power(a,n/2+1);
        t = (t + t * cur % M) % M;
        t = (t + cur) % M;
    }
    else
    {
        LL cur = power(a,n/2);
        t = (t + t * cur % M) % M;
    }
    return t;
}

int main()
{
    LL a,n;
    while(cin>>a>>n)
      cout<<sum(a,n)<<endl;
    return 0;
}

题目连接：http://poj.org/problem?id=3233

分析：矩阵求和S=（A+A^2+A^3+...+A^n）mod M.

矩阵快速幂+二分求和。

实现代码如下：

#include <cstdio>
#include <iostream>
using namespace std;
#define maxn 35
typedef struct
{
    int m[maxn][maxn];
}matrix;
matrix a,per;
int n,k,m;
void init()
{
    for(int i=0;i<n;i++)
      for(int j=0;j<n;j++)
      {
          scanf("%d",&a.m[i][j]);
          a.m[i][j]%=m;
          per.m[i][j]=(i==j);
      }
}
matrix add(matrix a,matrix b)
{
    matrix c;
    for(int i=0;i<n;i++)
      for(int j=0;j<n;j++)
        c.m[i][j]=(a.m[i][j]+b.m[i][j])%m;
    return c;
}
matrix multi(matrix a,matrix b)
{
    matrix c;
    for(int i=0;i<n;i++)
      for(int j=0;j<n;j++)
      {
          c.m[i][j]=0;
          for(int k=0;k<n;k++)
            c.m[i][j]+=a.m[i][k]*b.m[k][j];
        c.m[i][j]%=m;
      }
      return c;
}
matrix power(int k)
{
    matrix c,p,ans=per;
    p=a;
    while(k)
    {
        if(k&1)
        {
            ans=multi(ans,p);
            k--;
        }
        else
        {
            k/=2;
            p=multi(p,p);
        }
    }
    return ans;
}
matrix matrix_sum(int k)
{
    if(k==1) return a;
    matrix temp,b;
    temp=matrix_sum(k/2);
    if(k&1)
    {
        b=power(k/2+1);
        temp=add(temp,multi(temp,b));
        temp=add(temp,b);
    }
    else
    {
        b=power(k/2);
        temp=add(temp,multi(temp,b));
    }
    return temp;
}
int main()
{
    while(scanf("%d%d%d",&n,&k,&m)!=-1)
    {
        init();
        matrix ans=matrix_sum(k);
        for(int i=0;i<n;i++)
        {
            for(int j=0;j<n-1;j++)
              printf("%d ",ans.m[i][j]);
            printf("%d\n",ans.m[i][n-1]);
        }
    }
    return 0;
}