POJ 3233-Matrix Power Series 【快速幂】

Description

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 2
2 3

题意：

已知n，k，m，S = A + A2 + A3 + … + Ak. 求：S%m

思路1：

快速幂取模+二分优化。
当k为偶数时：
比方 k=6,那么 A+A^2+A^3+A^4+A^5+A^6 = A+A^2+A^3+ A^3*（A+A^2+A^3）
s(k) = s(k/2)+A^(n/2) * s(k/2) 即s(k) = (E+A^(n/2))*s(n/2) （E为单位矩阵）
当k为奇数时：
s(k)=s(k-1)+A^k , 那么k-1为偶数。能够依照上面的二分

ac代码1：

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <string>
#include <algorithm>
#include <cmath>
#include <functional>
#include <stack>
#include <queue>
#include <vector>
#include <map>
#include <set>
//#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll maxn = 81;
ll n, mod;
struct Mat {
    ll m[maxn][maxn];
};
Mat a, per;
void init()
{
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
        {
            cin >> a.m[i][j];
            a.m[i][j] %= mod;
            per.m[i][j] = (i == j);
        }
}

Mat mul(Mat A, Mat B)
{
    Mat ans;
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
        {
            ans.m[i][j] = 0;
            for (int k = 0; k < n; k++)
                ans.m[i][j] += (A.m[i][k] * B.m[k][j]);
            ans.m[i][j] %= mod;
        }
    return ans;
}
Mat power(ll k)
{
    Mat p = a, ans = per;
    while (k)
    {
        if (k & 1)
        {
            ans = mul(ans, p);
            k--;
        }
        else
        {
            k = k >> 1;
            p = mul(p, p);
        }
    }
    return ans;
}

Mat add(Mat a, Mat b)
{
    Mat 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]) % mod;
    return c;
}
Mat sum(ll k)
{
    if (k == 1)
        return a;
    Mat temp, b;
    temp = sum(k / 2);
    if (k & 1)
    {
        b = power(k / 2 + 1);
        temp = add(temp, mul(temp, b));
        temp = add(temp, b);
    }
    else
    {
        b = power(k / 2);
        temp = add(temp, mul(temp, b));
    }
    return temp;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    ll k;
    while (cin >> n >> k >> mod)
    {
        init();
        Mat ans = sum(k);
        for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j)
                printf("%lld%c", ans.m[i][j], (j == n - 1) ? '\n' : ' ');
    }
    return 0;
}

思路2：

推理发现：Fn = A + A*F（n-1）然后我们可以构造矩阵：

那么我们就可以用一个矩阵快速幂了。

ac代码2：

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <string>
#include <algorithm>
#include <cmath>
#include <functional>
#include <stack>
#include <queue>
#include <vector>
#include <map>
#include <set>
//#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

int n, k, mod;//n阶方阵，k次幂
struct er{
    int hang;
    int lie;
    ll m[66][66];
}a, g, h, result;

er multi(er x1, er x2)//x1 * x2
{
    er tmp;
    tmp.hang = x1.hang;
    tmp.lie = x2.lie;
    for (int i = 0; i < tmp.hang; ++i)
    {
        for (int j = 0; j < tmp.lie; ++j)
        {
            (tmp.m[i][j]) = 0;
            for (int k = 0; k < x1.lie; ++k)
                (tmp.m[i][j]) = ( (tmp.m[i][j]) + (x1.m[i][k] * x2.m[k][j]) ) % mod;
        }
    }
    return tmp;
}

er qmi(er x, int y)
{
    er ans;/*2n单位矩阵*/
    ans.hang = 2*n;
    ans.lie = 2*n;
    for (int i = 0; i < 2*n; ++i)
    {
        for (int j = 0; j < 2*n; ++j)
        {
            if(i==j)
                ans.m[i][j] = 1;
            else
                ans.m[i][j] = 0;
        }
    }
    while (y)
    {
        if (y&1)
            ans = multi(ans, x);
            //ans = ans * x;
        x = multi(x, x);
        //x = x * x;
        y = y >> 1;
    }
    return ans;
}

int main()
{
    scanf("%d %d %d", &n, &k, &mod);
    a.hang = n   ,  a.lie = n;
    h.hang = n   ,  h.lie = 2*n;
    g.hang = 2*n ,  g.lie = 2*n;
    for (int i = 0; i < n; ++i)
    {                                   //g矩阵为      (a，0)
        for (int j = 0; j < n; ++j)     //              (a，1)
        {
            scanf("%lld", &a.m[i][j]);  //h矩阵为      (a, 1)
            g.m[i][j] = a.m[i][j];
            g.m[i+n][j] = a.m[i][j];
            g.m[i][j+n] = 0;
            if(i == j)
                g.m[i+n][j+n] = 1;
            else
                g.m[i+n][j+n] = 0;

            h.m[i][j] = a.m[i][j];
            if(i == j)
                h.m[i][j+n] = 1;
            else
                h.m[i][j+n] = 0;
        }
    }
    result = multi( h , qmi(g, k-1) );
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            printf("%lld%c", result.m[i][j], (j == n - 1) ? '\n' : ' ');
    return 0;
}