POJ - 3233 ——Matrix Power Series (矩阵乘矩阵)

本文介绍了一种使用矩阵快速幂解决特定数学问题的方法:给定一个n×n矩阵A及正整数k,求矩阵A从1次方到k次方的累加和,并对其取模运算。文中提供了两种实现思路及对应的代码示例。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

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

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). 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 22 3

思路:因为A+A^2+A^3......A^k可化为F(n)=A+A*F(n-1),所以可以变为
{F(k)}={F(n-1)}*{A,A}={A}*{A,A}^k-1
{  1 }    {   1  }    {0,1}   {1}   {0,1}
或者用二分A+A^2+A^3......A^k=(A+A^2+A^3......A^k/2)*k/2*(A+A^2+A^3......A^k/2);
这题我用的是矩阵套矩阵的方法,但是不知道wa在哪里,怎么测都是对的。估计是我太菜;

祥神32msAC代码:

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <map>
#include <math.h>
#define LL long long
#define FOR(i) for(int i=1;i<=n;i++)
using namespace std;
int d[100][100];
int m;
struct matrix
{
    int n;
    int a[100][100];
    matrix(int nn,int flag):n(nn)
    {
        if(flag==0) FOR(i)FOR(j)a[i][j]=0;
        if(flag==1) FOR(i)FOR(j)a[i][j]=i==j;
        if(flag==2) FOR(i)FOR(j)a[i][j]=d[i][j];
    }
    matrix operator *(const matrix &b)const
    {
        matrix c(n,0);
        FOR(i)FOR(j)if(a[i][j])
            FOR(k)c.a[i][k]+=a[i][j]*b.a[j][k];
        FOR(i)FOR(j)if(c.a[i][j]>m)c.a[i][j]%=m;
        return c;
    }
    matrix powmat(int x)
    {
        matrix ans(n,0),tmp(n,2);
        for(int i=n/2+1;i<=n;i++)
            ans.a[i-n/2][i]=1;
        for(; x; tmp=tmp*tmp,x>>=1)
            if(x&1)ans=ans*tmp;
        return ans;
    }
};
int n;
int main()
{
    int k;
    while(~scanf("%d%d%d",&n,&k,&m))
    {
        int g;
        memset(d,0,sizeof(d));
        for(int i=1; i<=n; i++)
        {
            d[i][i]=1;
            for(int j=1; j<=n; j++)
            {
                scanf("%d",&g);
                if(g>m)g%=m;
                d[i+n][j]=d[i+n][j+n]=g;
            }
        }
        matrix w=matrix(2*n,2).powmat(k);
        for(int i=1; i<=n; i++)
        {
            for(int j=1; j<=n; j++)
                printf("%d ",w.a[i][j]%m);
            printf("\n");
        }
    }
}

自己的找不到错误的wa代码(过些时间在回来重新做一遍):
#include<algorithm>
#include<string.h>
#include<stdio.h>
#define LL long long
using namespace std;
int n,m;
struct node
{
    LL f[35][35];
};
void init(node &a,node &b)
{
    for(int i=1; i<=n; i++)
    {
        b.f[i][i]=1;
    }
    memset(a.f,0,sizeof(a.f));
}
void jzcj(node &c,node &a,node &b)
{
    for(int i=1; i<=n; i++)
    {
        for(int j=1; j<=n; j++)
        {
            for(int k=1; k<=n; k++)
            {
                c.f[i][j]+=a.f[i][k]*b.f[k][j];
                c.f[i][j]%=m;
            }
        }
    }
}
void see(node a)
{
    for(int i=1; i<=n; i++)
    {
        for(int j=1; j<=n; j++)
        {
            printf("%d ",a.f[i][j]);
        }
        printf("\n");
    }
}
node jjzm(node a,int k)
{
    node p[3][3],z[3][3],v[3][3],ans;
    p[1][1]=p[1][2]=a;
    init(p[2][1],p[2][2]);
    z[1][1]=z[2][2]=p[2][2];
    z[1][2]=z[2][1]=p[2][1];
    while(k)
    {
        if(k&1)
        {
            for(int i=1; i<=2; i++)
            {
                for(int j=1; j<=2; j++)
                {
                    memset(v[i][j].f,0,sizeof(v[i][j].f));
                    for(int k=1; k<=2; k++)
                    {
                        jzcj(v[i][j],p[i][k],z[k][j]);
                    }
                }
            }
            for(int i=1; i<=2; i++)
            {
                for(int j=1; j<=2; j++)
                {
                    z[i][j]=v[i][j];
                }
            }
        }
        k/=2;
        for(int i=1; i<=2; i++)
        {
            for(int j=1; j<=2; j++)
            {
                memset(v[i][j].f,0,sizeof(v[i][j].f));
                for(int k=1; k<=2; k++)
                {
                    jzcj(v[i][j],p[i][k],p[k][j]);
                }
            }
        }
        for(int i=1; i<=2; i++)
        {
            for(int j=1; j<=2; j++)
            {
                p[i][j]=v[i][j];
            }
        }
    }
    memset(ans.f,0,sizeof(ans.f));
//    see(z[1][1]);
    jzcj(ans,a,z[1][1]);
    for(int i=1; i<=n; i++)
    {
        for(int j=1; j<=n; j++)
        {
            ans.f[i][j]+=z[1][2].f[i][j];
            ans.f[i][j]%=m;
        }
    }
    return ans;
}
int main()
{
    int k;
    scanf("%d%d%d",&n,&k,&m);
    {
        node a,c;
        for(int i=1; i<=n; i++)
        {
            for(int j=1; j<=n; j++)
            {
                scanf("%I64d",&a.f[i][j]);
            }
        }
        c=jjzm(a,k-1);
        for(int i=1; i<=n; i++)
        {
            for(int j=1; j<=n; j++)
            {
                printf("%I64d ",c.f[i][j]%m);
            }
            printf("\n");
        }
    }
}



评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值