关灯游戏

例子:  

有一个5*6的灯泡构成的矩阵，灯的开关规则是这样：当改变某盏灯的，状态时，这盏灯的上下左右相邻的灯的状态也随之改变。例如：
0 1 1 0 1 0
1 0 0 1 1 1
0 0 1 0 0 1
1 0 0 1 0 1
0 1 1 1 0 0

当按下2行3列的开关时，状态变为：
0 1 0 0 1 0
1 1 1 0 1 1
0 0 0 0 0 1
1 0 0 1 0 1
0 1 1 1 0 0

游戏的目的是对于任意给定的亮灭初始态，通过一系列动作关闭所有的灯。
可以注意到的是：
1.矩阵的状态与按开关的顺序无关
2.如果某个开关按下了两次，那么就相当于取消了第一次的操作，也就是说没有开关需要按超过1次

现在问题是：对于给定的初始状态，求出需要按哪些开关来完成游戏

原题在这里
http://acm.pku.edu.cn/JudgeOnline/showproblem?problem_id=1222

把上面的矩阵看成一个m*n的向量X=(x1,x2,...,x(m*n))
对于位置k上的开关，它将变化最多5个位置的开关，对应一个向量
C(k)=(0,0,...,1,0,....,1,...,0)
其中开关状态改变的位置为1，开关状态不改变的位置为0
对于初始向量X=(x1,x2,...,x(m*n)),使用了开关C(k)后，状态会变成
X+C(k) (mod 2)

所以对初始向量X,我们需要选择一系列的k1,k2,...,ks使得
X+C(k1)+C(k2)+....+C(ks) (mod 2)=O=(0,0,0,...,0)
我们可以同样构造一个0,1向量Y,使得，如果位置k出现在k1,k2,...ks中，那么Y
在位置k的值是1，不然是0，这样，我们就可以将上面公式写成矩阵形式
X+Y*C (mod 2)=O
其中C=(C(1)' C(2)' .... C(m*n)')'
也就是C是由这m*n个行向量构成的矩阵，第k行就是向量C(k)
在二阶域上，加和减是相同的，也就是上面的方程等价于
Y*C (mod 2)=X
其中C,X已知，求Y.
由于(mod 2)运算是一个域 （关于乘除加减封闭，加减是mod 2加减，还满足结合率，交换率）
所以我们可以直接在二阶域上用高斯消元法求解（注意加减是mod 2的，对应计算机上的异或运算）
其中，如果C可逆，解是唯一的，如果C不可逆，解可能不存在，也可能不唯一。而如果解不唯一时，

这时解德数目同C的秩有关系，比较有意思的是，如果我们记f(n,2x)是n次第二类切皮雪夫多项式，

那么gcd(f(m,x)(mod 2),f(n,1-x)(mod 2))的次数就是C的自由度（也就是矩阵的阶数减去秩）。

具体证明可以在百度贴吧中内容。

而这里C是一个分块三对角阵，我们可以使用追逐法求解,时间复杂度可以只有O(n^3).比较有意思的是这里的分析

同另外一个扫雷的算法问题几乎完全相同，区别是那里在普通的数域，而这里在二阶有限域上。

下面贴出的程序的时间复杂度为O(n^4),这是因为其中同稀疏矩阵的乘法没有优化.

源代码和预编译的代码可以在
     http://download.csdn.net/source/3275909
下载
/*The code to turn off all bulbs in a matrix of bulbs
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

typedef char ** MATRIX;
typedef char * VECTOR;
typedef char *const* CONST_MATRIX;
typedef const char * CONST_VECTOR;

MATRIX matrix_alloc2(int n, int m){
 MATRIX x= (MATRIX)malloc(sizeof(char *)*n+sizeof(char)*n*m);
 int i;
 x[0]=(char *)(x+n);
 for(i=1;i<n;i++)x[i]=x[i-1]+m;
 return x;
}

int maxtrix_count_one(MATRIX x,int n, int m){
 int r=0;
 int i,j;
 for(i=0;i<n;i++)for(j=0;j<n;j++)if(x[i][j]!=0)r++;
 return r;
}

MATRIX matrix_alloc(int n){
    MATRIX x = (MATRIX)malloc(sizeof(char *)*n+sizeof(char)*n*n);
    int i;
    x[0]=(char *)(x+n);
    for(i=1;i<n;i++)x[i]=x[i-1]+n;
    return x;
}

void matrix_free(MATRIX x){
    free(x);
}

VECTOR vector_alloc(int n){
    VECTOR x = (VECTOR)malloc(sizeof(char)*n);
    return x;
}

void vector_free(VECTOR x){
    free(x);
}

void matrix_sum(MATRIX A, CONST_MATRIX B, int n){
    int i,j;
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            A[i][j] ^= B[i][j];
        }
    }
}

void vector_sum(VECTOR a, CONST_VECTOR b, int n){
    int i;
    for(i=0;i<n;i++)a[i]^=b[i];
}

void matrix_mul(MATRIX out, CONST_MATRIX in1, CONST_MATRIX in2, int n){
    int i,j,k;
    for(i=0;i<n;i++)for(j=0;j<n;j++){
        int sum=0;
        for(k=0;k<n;k++){
            if(in1[i][k])
                sum^=in2[k][j];
        }
        out[i][j]=sum;
    }
}

void matrix_mul_H(MATRIX out, CONST_MATRIX in, int n){
    int i,j;
    for(i=0;i<n;i++)for(j=0;j<n;j++){
        int sum=in[i][j];
        if(i>0)sum^=in[i-1][j];
        if(i<n-1)sum^=in[i+1][j];
        out[i][j]=sum;
    }
}

void matrix_mul_vector(VECTOR out, CONST_MATRIX m, CONST_VECTOR in, int n){
    int i,j;
    for(i=0;i<n;i++){
        int sum=0;
        for(j=0;j<n;j++){
            sum^=m[i][j]&in[j];
        }
        out[i]=sum;
    }
}

void H_mul_vector(VECTOR out, CONST_VECTOR in, int n){
    int i;
    for(i=0;i<n;i++){
        int sum=in[i];
        if(i>0)sum^=in[i-1];
        if(i<n-1)sum^=in[i+1];
        out[i]=sum;
    }
}

void vector_init_const(VECTOR v, char c, int n){
    int i;
    for(i=0;i<n;i++)v[i]=c;
}

void matrix_init_O(MATRIX o, int n){
    int i,j;
    for(i=0;i<n;i++)for(j=0;j<n