矩阵和线性方程组

总结：满足下面这句话的条件，都可以考虑矩阵乘法+快速幂——把一个向量v变成另一个向量v'，并且v'的每一个分量都是v的各个分量的线性组合。

 

1、UVa 10870, Recurrences

题意：给出递推关系f(n) = a1*f(n-1) + a2*f(n-2)....ad*f(n-d)，给出d, n, m, a1, a2, ...ad, f(1), f(2)...f(n-d)，求f(n) % m。

解法：矩阵乘法+快速幂。这是我第一次做矩阵乘法的题，由于是专题训练，一上来就构造成f(n) = G * F(n-1)，只需要G = [a1, a2,...ad]，F(n-1) = [f(n-1), f(n-2),...f(n-d)]T即可。

　　　但是，这样构造显然是不对的，我们要构造出的灯式形势应该是这样的(矩阵)F(n) = (d*d的矩阵)G * (矩阵)F(n-1)，只有这样的形式(G为d*d)，才能使用矩阵乘法快速幂。

　　　所以，F(n) = [f(n-d), f(n-d+1),...f(n-1)]T，

　　　G = [0   1　　　　　　　　]

    　　　   [    0    1 　　　　     ]

　　　　　 [...　　　　　...　　　 ]

　　　　　 [                     0     1]　　　　　　　　　　　　　　　　　　　　　　

　　　　　 [ad a(d-1)....    a2  a1]　(d*d的矩阵)

　　　构造出来，这道题便得以解决。

tag:math, matrix

/*
 * Author:  Plumrain
 * Created Time:  2013-09-10 21:16
 * File Name: math-UVa-10870.cpp
 */
#include<iostream>
#include<cstdio>
#include<cstring>

using namespace std;

#define CLR(x) memset(x, 0, sizeof(x))
#define out(x) cout<<#x<<":"<<(x)<<endl
#define tst(a) cout<<#a<<endl

typedef int matrix[20][20];

const int N = 15;
int mod;
int a[N+5], f[N+5];

void mtx_init (matrix& A, int n)
{
    CLR (A);
    for (int i = 0; i < (n-1); ++ i)
        A[i][i+1] = 1;
    for (int i = 0; i < n; ++ i)
        A[n-1][i] = a[n-i];
}

void mtx_mul (matrix& A, matrix B, int n)
{
    matrix ret;
    CLR (ret);
    for (int i = 0; i < n; ++ i)
        for (int k = 0; k < n; ++ k)
            for (int j = 0; j < n; ++ j)
                ret[i][k] = (((A[i][j] * B[j][k]) % mod) + ret[i][k]) % mod;

    for (int i = 0; i < n; ++ i)
        for (int j = 0; j < n; ++ j)
            A[i][j] = ret[i][j];
}
void mtx_pow (matrix& A, int num, int n)
{
    matrix ret;
    CLR (ret);
    for (int i = 0; i < n; ++ i)
        ret[i][i] = 1; 
    while (num > 0){
        if (num & 1) 
            mtx_mul (ret, A, n);
        num >>= 1;
        mtx_mul (A, A, n);
    }

    for (int i = 0; i < n; ++ i)
        for (int j = 0; j < n; ++ j)
            A[i][j] = ret[i][j];
}

int gao (int m, int n)
{
    matrix A;
    mtx_init (A, n);
    
    mtx_pow (A, m-n, n); 
    int ret = 0;
    for (int i = 0; i < n; ++ i)
        ret = (((A[n-1][i] * f[i+1]) % mod) + ret) % mod;
    if (ret < 0) ret += mod;
    return ret;
}

int main()
{
    int n, m; 
    while (scanf ("%d%d%d", &n, &m, &mod) != EOF && n){
        for (int i = 1; i <= n; ++ i){
            scanf ("%d", &a[i]);
            a[i] %= mod;
        }
        for (int i = 1; i <= n; ++ i){
            scanf ("%d", &f[i]);
            f[i] %= mod;
        }

        if (m <= n) printf ("%d\n", f[m]);
        else printf ("%d\n", gao (m, n));
    }
    return 0;
}

 

2、NEERC 2006, LA 3704, Cellular Automaton

题意：一个环上有n个点，标号依次为0,1...n-1。每次操作后，每个点的值都将变为，与它距离小于等于d的所有点在操作之前的值之和mod m的值。给定n, m, d, k和n个格子各自的值，求经过k次操作之后每个格子的值。

解法：

tag:math, circulant matrix,

/*
 * Author:  Plumrain
 * Created Time:  2013-09-13 12:59
 * File Name: math-NEERC2006-LA3704.cpp
 */
#include<iostream>
#include<cstdio>
#include<cstring>

using namespace std;

#define CLR(x) memset(x, 0, sizeof(x))

const int N = 500;

typedef long long int64;
typedef int64 matrix[N+5];

int n, mod, d, k;
int64 a[N+5];
int64 A[N+5];

void mtx_mul(matrix& A, matrix B)
{
    matrix ret;
    CLR (ret);
    for (int i = 0; i < n; ++ i)
        for (int j = 0, pos = i; j < n; ++ j, pos=(pos-1+n)%n)
            ret[i] = ((A[j] * B[pos]) % mod + ret[i]) % mod;

    for (int i = 0; i < n; ++ i)
        A[i] = ret[i];
}

void mtx_pow(matrix& A, int num)
{
    if (num == 1) return;

    matrix ret;
    for (int i = 0; i < n; ++ i)
        ret[i] = A[i];
    -- num;

    while (num > 0){
        if (num & 1) mtx_mul (ret, A);
        num >>= 1;
        mtx_mul (A, A);
    }

    for (int i = 0; i < n; ++ i)
        A[i] = ret[i];
}

void gao()
{
    CLR (A);
    int pos1 = 0, pos2 = 0;
    int cnt = 0;
    while (cnt <= d){
        A[pos1] = 1; A[pos2] = 1;

        ++ cnt;
        pos1 = (pos1+1) % n;
        pos2 = (pos2-1+n) % n;
    }

    mtx_pow(A, k);

    for (int i = 0; i < n; ++ i){
        int64 x = 0;
        for (int j = 0; j < n; ++ j)
            x = ((A[(j-i+n)%n]*a[j]) % mod + x) % mod;

        printf ("%lld", x);
        if (i != n-1) printf (" ");
    }
    printf ("\n");
}

int main()
{
    while (scanf ("%d", &n) != EOF){    
        scanf ("%d%d%d", &mod, &d, &k);
        for (int i = 0; i < n; ++ i)
            scanf ("%lld", &a[i]);

        gao ();
    }
    return 0;
}

 

3、POJ 3070 Fibonacci

题意：对于Fibonacci数列，“0, 1, 1, 2, 3, 5, 8, 13, 21, 34, …”，(F0 = 0)可以用矩阵的方法求得Fn。给n，求Fn。

解法：本来，如果直接说给n求Fn，这道题还是很有价值的。但是题目不仅仅是直接说了用矩阵求，甚至直接将给出如下图片：

这道题就直接按照这个图片写就好了.....

Ps：用矩阵求类似问题的构造方法，还有就是刘汝佳书上写的那种构造方法，那个是通用的。即本篇文章第一题。

tag:math, matrix, good

/*
 * Author:  Plumrain
 * Created Time:  2013-10-14 15:42
 * File Name: math-POJ-3070.cpp
 */
#include<iostream>
#include<cstdio>
#include<cstring>

using namespace std;

#define CLR(x) memset(x, 0, sizeof(x))
const int mod = 10000;
typedef int matrix[10][10];

void mtx_mul(matrix& A, matrix B)
{
    matrix ret; CLR (ret);
    for (int i = 0; i < 2; ++ i)
        for (int j = 0; j < 2; ++ j){
            ret[i][j] = 0;
            for (int k = 0; k < 2; ++ k)
                ret[i][j] = (A[i][k] * B[k][j] + ret[i][j]) % mod;
        }

    for (int i = 0; i < 2; ++ i)
        for (int j = 0; j < 2; ++ j)
            A[i][j] = ret[i][j];
}

void mtx_pow(matrix& A, int n)
{
    -- n;
    matrix ret; CLR (ret);
    for (int i = 0; i < 2; ++ i)
        for (int j = 0; j < 2; ++ j)
            ret[i][j] = A[i][j];
    while (n > 0){
        if (n & 1) mtx_mul(ret, A);
        n >>= 1;
        mtx_mul(A, A);
    }
    for (int i = 0; i < 2; ++ i)
        for (int j = 0; j < 2; ++ j)
            A[i][j] = ret[i][j];
}

int main()
{
    int n;
    while (scanf ("%d", &n) != EOF && n != -1){
        matrix A; CLR (A);
        A[0][0] = 1; A[0][1] = 1; A[1][0] = 1;
        if (n < 4){
            if (!n) printf ("0\n");
            if (n == 1 || n == 2) printf ("1\n");
            if (n == 3) printf ("2\n");
            continue;
        }
        mtx_pow(A, n);
        printf ("%d\n", A[0][1] % mod);     
    }
    return 0;
}