数论之矩阵快速幂和二分法 poj3233

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 10⁹) and m (m < 10⁴). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order.

Output

Output the elements of S modulo m in the same way as A is given.

Sample Input

2 2 4
0 1
1 1

Sample Output

1 2
2 3

求和时利用二分思想：当n为奇数时，f(n)=f(n-1)+A^n,当n时偶数时，f(n)=f(n/2)+f(n/2)*A^(n/2);

矩阵相乘时注意优化剪枝，求矩阵的N次方时也用到了二分思想；

代码如下：

#include<stdio.h>
#include<iostream>
#include<string.h>
using namespace std;
int n,m,k;
struct dog
{
    int n[35][35];
};
dog m1;
dog mul(dog a,dog b)
{
    dog ret;
    int i,j,z;
    memset(ret.n,0,sizeof(ret.n));
    for(i=0;i<n;++i)
    {
        for(z=0;z<n;++z)
        if(a.n[i][z])
        for(j=0;j<n;++j)
        {
            ret.n[i][j]+=a.n[i][z]*b.n[z][j];
            if(ret.n[i][j]>=m)ret.n[i][j]%=m;
        }
    }
    return ret;
}
/*dog exp(dog a, int k)
{
    if(k==1)return a;
    dog e;
    memset(e.n,0,sizeof(e.n));
    for(int i=0;i<n;++i){e.n[i][i]=1;}
    if(k==0)return e;
    while(k)
    {
        if(k&1)e=mul(a,e);
        a=mul(a,a);
        k>>=1;
    }
    return e;
}*/
dog exp(dog a,int l)
{
    if(l==1)
    return a;
    dog e;
    memset(e.n,0,sizeof(e.n));
    for(int i=0;i<n;i++)
    e.n[i][i]=1;
    if(l==0)
    return e;
    while(l)
    {
        if(l&1)
        e=mul(a,e);
        a=mul(a,a);
        l>>=1;
    }
    return e;
}
dog add(dog a,dog b)
{
    int i,j;
    dog c;
    memset(c.n,0,sizeof(c.n));
    for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    {
        c.n[i][j]=a.n[i][j]+b.n[i][j];
        if(c.n[i][j]>=m)
        c.n[i][j]%=m;
    }
    return c;
}
dog sum(int t)
{
    if(t==1)
    return m1;
    if(t&1)
    return add(sum(t-1),exp(m1,t));
    else
    {
        dog s=sum(t>>1);
        return add(s,mul(s,exp(m1,t>>1)));
    }

}
int main()
{
    dog sum1;
    memset(sum1.n,0,sizeof(sum1.n));
    while(~scanf("%d%d%d",&n,&k,&m))
    {
    for(int i=0;i<n;i++)
    for(int j=0;j<n;j++)
    {
        scanf("%d",&m1.n[i][j]);
        if(m1.n[i][j]>=m)
        m1.n[i][j]%=m;
    }
    sum1=sum(k);
    for(int w=0;w<n;w++)
    {
        printf("%d",sum1.n[w][0]);
        for(int j=1;j<n;j++)
        printf(" %d",sum1.n[w][j]);
        printf("\n");
    }
    }

    return 0;
}