有理数域上lamada矩阵(方阵)化标准形的实现

前言：首先介绍一下lamada矩阵，其为高等代数或线性代数的内容。其中将λ-矩阵化成标准形在这门课中占据着举足轻重的地位。lamada矩阵即λ-矩阵，亦称多项式矩阵，是以多项式为元素的矩阵。而今天要研究的就是在有理数域上的多项式组成的矩阵。而数字矩阵是λ-矩阵的一种特殊情形。下面介绍λ-矩阵的标准形，即除主对角线上的元素外其余元素全为0，且主对角线上的前一个元素总是能整除后一个元素。

据高代代数的知识，矩阵化标准形主要是通过几个初等变换，即换法变换，倍法变换以及消法变换。而λ-矩阵化标准形除需要以上知识外还需要一个引理：

若一个λ-矩阵的左上角元素a11不为0，记为A，且A中至少有一个不能被a11整除，那么就一定能够通过初等变换找到与之相似的另一个λ-矩阵B,其中B的左上角元素a‘11也不为0，但其次数比a11低。

那么，最头疼的问题来了，该怎么表示一个λ-矩阵呢？wo苦思冥想，翻来覆去，痛定思痛，奋起反抗，前后尝试了链表加链表，链表加数组，链表加链表加链表的形式存储，但都因为某些莫名奇妙无法正常赋值的原因失败了。所以我决定开摆，二维数组加二维数组，通过映射连接，除了不太好赋值之外，其他都还好。举个栗子。比如下面这个λ-矩阵。

就可以这么表示：

其中n是λ-矩阵的阶数，DEG是多项式最大的次数。定义如下：

一路艰难，咱到这算是能够表示一个λ-矩阵了。接下来我们着重解决化λ-矩阵为标准形的准备工作，包括：

三个初等初等变换的实现

以引理为理论基础通过初等变换将λ-矩阵左上角元素化成能整除所有元素

当然还有一些对于多项式的一些基础操作咱得先做好准备。

例如，多项式的更新：

//改变多项式，输入原多项式，新多项式 
void modify(int *prime,int *assig)
{for(int i=0;i<DEG;i++){prime[i] = assig[i];}
}

再如，多项式的次数：

//获取数组长度，即多项式次数 
int strlenNum(int *p)
{int k=DEG-1,res;while(p[k]==0 && k>=0){k--;}res = k + 1;return res;
}

再如，多项式加减法：

//多项式加减法，输入两个多项式，q加(减)p, 结果保存至传入的指针变量，以及flag1表示加，0表示减 
void addsub(int *q,int *p,int *result,int flag)
{int strq = strlenNum(q),strp = strlenNum(p);int Numax = max(strq,strp);for(int i=0;i<DEG;i++) result[i]=0;if(flag){for(int i=0;i<Numax;i++){result[i] = q[i] + p[i];}}else{for(int i=0;i<Numax;i++){result[i] = q[i] - p[i];}}
}

再如，多项式乘法：

//多项式乘法，输入两个多项式，将相乘结果保存至传入的指针变量
void polymulti(int *q,int *p,int *result)
{int strq = strlenNum(q),strp = strlenNum(p);for(int i=0;i<DEG;i++) result[i]=0;int i,j;for(i=0;i<strq;i++){for(j=0;j<strp;j++){result[i+j] += q[i]*p[j];}}
}

又如多项式的除法：

//多项式除法，输入两个多项式m除n，输出q(x),r(x)，能整除输出1，不能输出0 
int division(int *s,int *t,int *q,int *r)
{int N[DEG];//获取数组长度，即多项式系数 int strm = strlenNum(s),strn = strlenNum(t);if(strm > 0){//获取对n多项式作变换的倍数 int flag1 = strn - strm + 1,flag2 = s[strm-1];int temp[DEG] = {flag1*flag2};polymulti(t,temp,N);//准备完毕，开始相除//q(x),r(x)初始赋值for(int i=0;i<DEG;i++){q[i] = 0;r[i] = 0;}//定义除法产生的临时多项式 int mm[DEG],nn[DEG],bb[DEG];//初始化 modify(r,N);int a = strn - strm,b = strn - 1;while(strlenNum(r) >= strm && a>=0 && b>=0){q[a--] = r[b--] / s[strm-1];for(int i=0;i<DEG;i++) bb[i] = 0;bb[a+1] = q[a+1]; polymulti(s,bb,mm);//相乘 addsub(r,mm,nn,0);//相减 modify(r,nn);}if(!strlenNum(r)) return 1;else return 0;}else{exit(-1);}
}

对于多项式的除法，因为C语言对于非整形数的处理相当麻烦，而且，将多项式放大非0常数倍并不会改变两个多项式之间的整除关系。这里我们采用的是将被除数放大flag1*flag2倍的方式以达到所有运算都是整数运算的目的(其中呢flag1是被除多项式次数减除多项式的次数加一，flag2是除多项式系数最高次数系数)。

至此呢，我们的一些对于多项式的基础操作完成了，下面开始三个初等变换的操作：

首先是换法变换，就是将第i行(列)和第j行(列)互换，很简单：

//换法变换，输入lamada矩阵，需要换的行(列)i，j，以及判断变量1表示换行，0表示换列 
void change(int (*lamada)[n],int i,int j,int flag)
{int temp;i = i-1;j = j-1;if (flag) {for(int k=0;k<n;k++){temp = lamada[i][k];lamada[i][k] = lamada[j][k];lamada[j][k] = temp;}} else {for(int k=0;k<n;k++){temp = lamada[k][i];lamada[k][i] = lamada[k][j];lamada[k][j] = temp;}}
}

接着是倍法变换，就是将第k行(列)乘一个非0的多项式，也不难，主要用到多项式的乘法：

//倍法变换，输入lamada矩阵，原系数矩阵，倍数为一个多项式，行(列)， 以及判断变量1表示换行，0表示换列
void multiple(int (*lamada)[n],int (*coef)[DEG],int *times,int k,int flag)
{int res[DEG];k = k - 1;if (flag) {for(int i=0;i<n;i++){polymulti(coef[lamada[k][i]],times,res);modify(coef[lamada[k][i]],res);}} else {for(int i=0;i<n;i++){polymulti(coef[lamada[i][k]],times,res);modify(coef[lamada[i][k]],res);}}
}

最后是消法变换，也是和化标准形最接近的一个变换方式，主要用到多项式乘法和加减法：

//消法变换，输入lamada矩阵，原系数矩阵，第i行(列)减去第j行(列)的times倍，以及判断变量1表示换行，0表示换列 
void remove(int (*lamada)[n],int (*coef)[DEG],int i,int j,int *times,int flag)
{int temp1[DEG],temp2[DEG];i = i - 1;j = j - 1;if (flag) {for(int k=0;k<n;k++){polymulti(coef[lamada[j][k]],times,temp1);addsub(coef[lamada[i][k]],temp1,temp2,0);modify(coef[lamada[i][k]],temp2);}} else {for(int k=0;k<n;k++){polymulti(coef[lamada[k][j]],times,temp1);addsub(coef[lamada[k][i]],temp1,temp2,0);modify(coef[lamada[k][i]],temp2);}}
}

至此呢，第一个准备工作就完成了，是不是感觉很简单呢。如果是那就分享转发一下吧，如果感觉有点懵咱就收藏一下，以后多看几遍。ok打起精神，我们继续前进。

第二个任务是啥来着，回去看看？不需要，我告诉你，第二个任务是将λ-矩阵左上角元素化成能整除所有元素。用到那个引理，这里不太好解释，直接上代码吧。

//对于第一行(列)未能被整除的情况降次 输入1表示第一行有不能被整除的，0表示第一列有不能被整除的 
//index表示所要处理的第几个左上角元素 以0开始 k是某元素列(行)下标 
void rl_decline(int (*lamada)[n],int (*coef)[DEG],int k,int flag,int index)
{int s[DEG],t[DEG],Q[DEG],R[DEG];modify(s,coef[lamada[index][index]]);if(flag)modify(t,coef[lamada[index][k]]);elsemodify(t,coef[lamada[k][index]]);division(s,t,Q,R);int strs = strlenNum(s),strt = strlenNum(t);int flag1 = strt - strs + 1,flag2 = s[strs-1];int times[DEG] = {flag1*flag2};if(flag){for(int i=index;i<n;i++){modify(t,coef[lamada[i][k]]);polymulti(t,times,coef[lamada[i][k]]);}remove(lamada,coef,k+1,index+1,Q,0);change(lamada,k+1,index+1,0);}else{for(int i=index;i<n;i++){modify(t,coef[lamada[k][i]]);polymulti(t,times,coef[lamada[k][i]]);}remove(lamada,coef,k+1,index+1,Q,1);change(lamada,k+1,index+1,1);}
}
//对于非第一行(列)未能被整除的降次，x，y都为数组元素下标 
//index表示所要处理的第几个左上角元素 以0开始 
void all_decline(int (*lamada)[n],int (*coef)[DEG],int x,int y,int index)
{int temp1[DEG],R[DEG],deone[DEG]={-1};int flag = division(coef[lamada[index][index]],coef[lamada[x][0]],temp1,R);if(!flag) return;remove(lamada,coef,x+1,index+1,temp1,1);remove(lamada,coef,index+1,x+1,deone,1);rl_decline(lamada,coef,y,1,index);
}
//初等变换使得左上角元素能整除所有元素，本质就是不断降低次数 
void ele_oper(int (*lamada)[n],int (*coef)[DEG],int index)
{int i,j,r_index=index,l_index=index,min = strlenNum(coef[lamada[index][index]]);int prem[DEG],Q[DEG],R[DEG];modify(prem,coef[lamada[index][index]]);for(i=index;i<n;i++){for(j=index;j<n;j++){if(strlenNum(coef[lamada[i][j]]) < min){min = strlenNum(coef[lamada[i][j]]);r_index=i;l_index=j;}}}change(lamada,index+1,r_index+1,1);change(lamada,index+1,l_index+1,0);for(i=index;i<n;i++){if(!division(prem,coef[lamada[i][index]],Q,R) || !division(prem,coef[lamada[index][i]],Q,R)){break;}}if(i<n && !division(prem,coef[lamada[index][i]],Q,R)){rl_decline(lamada,coef,i,1,index);}else if(i<n && !division(prem,coef[lamada[i][index]],Q,R)){rl_decline(lamada,coef,i,0,index);}else if(i==n){for(i=index+1;i<n;i++){for(j=index+1;j<n;j++){if(!division(prem,coef[lamada[i][j]],Q,R))r_index = i;l_index = j;}}all_decline(lamada,coef,r_index,l_index,index);}
}
//将左上角元素变成能整除所有元素 
void middle_temp(int (*lamada)[n],int (*coef)[DEG],int index)
{int len = strlenNum(coef[lamada[index][index]]);int flag = all_division(lamada,coef,index);if(len==1 || flag) return;else {ele_oper(lamada,coef,index);middle_temp(lamada,coef,index);}
}

大概思路就是第一行(列)上有不能被整除的元素，就能够通过函数rl_decline()降低左上角元素的次数，而当存在非第一行(列)的元素且第一行(列)上所有元素都能被整除时，可通过函数all_decline()将其化成第一行有不能被整除的元素，那么就可以通过函数rl_decline()降低左上角元素的次数。这样反反复复总能化成左上角元素能整除所有元素。

ok，虽说我不能解释这个第二个任务具体是怎么实现的，但是聪明的你一定能从我这屎一样的代码中看出来原理的对吧！

到这里呢，我们的准备工作就全部完成。下面才真正的进入到化标准形的过程。

首先呢，我们能将一个λ-矩阵化成左上角元素a11能整除所有元素的矩阵，那么我们通过“剥洋葱”的方法，通过消法变换，因为a11能整除所有的元素，那么就一定能找到一个多项式q使得第i行的第一个元素等于a11*q，那么就可以通过消法变换让第i行减去第一行的q倍，这样第i行的第一个元素就能化为0了，同样的对于每一行(大于1)的第一个元素都能通过这种方法化为0，那么再所有行(大于1)的第一个元素化为0之后，每一列的第一个元素，也能找到一个多项式p使得其等于a11*p，这样，同样的使用消法变换，每一列减第一列的那个特定的p倍就也能将每一列的第一个元素化为0，至此，最外层的多项式就全部华为0了。

接着可以发现，里面存在一个阶数减一的矩阵，同样重复以上所有步骤，找到左上角能整除所有的元素，再将最外层元素化为0，这样不断迭代，就可以得到一个λ-矩阵的标准形了。

下面是上述过程的代码实现：

//开始真正的标准形初等变换
void rel_eleoper(int (*lamada)[n],int (*coef)[DEG],int index)
{int strlen,flag1,flag2,Q[DEG],R[DEG],prem[DEG];modify(prem,coef[lamada[index][index]]);int zeros[DEG] = {0},strp = strlenNum(prem);for(int i=index+1;i<n;i++){strlen = strlenNum(coef[lamada[i][index]]);flag1 = strlen - strp + 1;flag2 = prem[strp-1];int temp[DEG] = {flag1*flag2};division(prem,coef[lamada[i][index]],Q,R);multiple(lamada,coef,temp,i+1,1);remove(lamada,coef,i+1,index+1,Q,1);}for(int i=index+1;i<n;i++){modify(coef[lamada[index][i]],zeros);}
}

最后呢给出整个代码，以及结果图(结果是开头说到的矩阵的标准形)：

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
//这里多项式的定义最大次数 
#define DEG 10
//这里定义lamada矩阵阶数 
#define n 3
//取两数最大
int max(int a,int b)
{if(a>b) return a;else return b;
}
//获取数组长度，即多项式次数 
int strlenNum(int *p)
{int k=DEG-1,res;while(p[k]==0 && k>=0){k--;}res = k + 1;return res;
}
//改变多项式，输入原多项式，新多项式 
void modify(int *prime,int *assig)
{for(int i=0;i<DEG;i++){prime[i] = assig[i];}
}
//换法变换，输入lamada矩阵，需要换的行(列)i，j，以及判断变量1表示换行，0表示换列 
void change(int (*lamada)[n],int i,int j,int flag)
{int temp;i = i-1;j = j-1;if (flag) {for(int k=0;k<n;k++){temp = lamada[i][k];lamada[i][k] = lamada[j][k];lamada[j][k] = temp;}} else {for(int k=0;k<n;k++){temp = lamada[k][i];lamada[k][i] = lamada[k][j];lamada[k][j] = temp;}}
}
//多项式乘法，输入两个多项式，将相乘结果保存至传入的指针变量
void polymulti(int *q,int *p,int *result)
{int strq = strlenNum(q),strp = strlenNum(p);for(int i=0;i<DEG;i++) result[i]=0;int i,j;for(i=0;i<strq;i++){for(j=0;j<strp;j++){result[i+j] += q[i]*p[j];}}
}
//倍法变换，输入lamada矩阵，原系数矩阵，倍数为一个多项式，行(列)， 以及判断变量1表示换行，0表示换列
void multiple(int (*lamada)[n],int (*coef)[DEG],int *times,int k,int flag)
{int res[DEG];k = k - 1;if (flag) {for(int i=0;i<n;i++){polymulti(coef[lamada[k][i]],times,res);modify(coef[lamada[k][i]],res);}} else {for(int i=0;i<n;i++){polymulti(coef[lamada[i][k]],times,res);modify(coef[lamada[i][k]],res);}}
}
//多项式加减法，输入两个多项式，q加(减)p, 结果保存至传入的指针变量，以及flag1表示加，0表示减 
void addsub(int *q,int *p,int *result,int flag)
{int strq = strlenNum(q),strp = strlenNum(p);int Numax = max(strq,strp);for(int i=0;i<DEG;i++) result[i]=0;if(flag){for(int i=0;i<Numax;i++){result[i] = q[i] + p[i];}}else{for(int i=0;i<Numax;i++){result[i] = q[i] - p[i];}}
}//消法变换，输入lamada矩阵，原系数矩阵，第i行(列)减去第j行(列)的times倍，以及判断变量1表示换行，0表示换列 
void remove(int (*lamada)[n],int (*coef)[DEG],int i,int j,int *times,int flag)
{int temp1[DEG],temp2[DEG];i = i - 1;j = j - 1;if (flag) {for(int k=0;k<n;k++){polymulti(coef[lamada[j][k]],times,temp1);addsub(coef[lamada[i][k]],temp1,temp2,0);modify(coef[lamada[i][k]],temp2);}} else {for(int k=0;k<n;k++){polymulti(coef[lamada[k][j]],times,temp1);addsub(coef[lamada[k][i]],temp1,temp2,0);modify(coef[lamada[k][i]],temp2);}}
}
//多项式除法，输入两个多项式m除n，输出q(x),r(x)，能整除输出1，不能输出0 
int division(int *s,int *t,int *q,int *r)
{int N[DEG];//获取数组长度，即多项式系数 int strm = strlenNum(s),strn = strlenNum(t);if(strm > 0){//获取对n多项式作变换的倍数 int flag1 = strn - strm + 1,flag2 = s[strm-1];int temp[DEG] = {flag1*flag2};polymulti(t,temp,N);//准备完毕，开始相除//q(x),r(x)初始赋值for(int i=0;i<DEG;i++){q[i] = 0;r[i] = 0;}//定义除法产生的临时多项式 int mm[DEG],nn[DEG],bb[DEG];//初始化 modify(r,N);int a = strn - strm,b = strn - 1;while(strlenNum(r) >= strm && a>=0 && b>=0){q[a--] = r[b--] / s[strm-1];for(int i=0;i<DEG;i++) bb[i] = 0;bb[a+1] = q[a+1]; polymulti(s,bb,mm);//相乘 addsub(r,mm,nn,0);//相减 modify(r,nn);}if(!strlenNum(r)) return 1;else return 0;}else{exit(-1);}
}
//至此，准备工作全部完毕，下面开始对lamada化标准形 
//对于第一行(列)未能被整除的情况降次 输入1表示第一行有不能被整除的，0表示第一列有不能被整除的 
//index表示所要处理的第几个左上角元素 以0开始 k是某元素列(行)下标 
void rl_decline(int (*lamada)[n],int (*coef)[DEG],int k,int flag,int index)
{int s[DEG],t[DEG],Q[DEG],R[DEG];modify(s,coef[lamada[index][index]]);if(flag)modify(t,coef[lamada[index][k]]);elsemodify(t,coef[lamada[k][index]]);division(s,t,Q,R);int strs = strlenNum(s),strt = strlenNum(t);int flag1 = strt - strs + 1,flag2 = s[strs-1];int times[DEG] = {flag1*flag2};if(flag){for(int i=index;i<n;i++){modify(t,coef[lamada[i][k]]);polymulti(t,times,coef[lamada[i][k]]);}remove(lamada,coef,k+1,index+1,Q,0);change(lamada,k+1,index+1,0);}else{for(int i=index;i<n;i++){modify(t,coef[lamada[k][i]]);polymulti(t,times,coef[lamada[k][i]]);}remove(lamada,coef,k+1,index+1,Q,1);change(lamada,k+1,index+1,1);}
}
//对于非第一行(列)未能被整除的降次，x，y都为数组元素下标 
//index表示所要处理的第几个左上角元素 以0开始 
void all_decline(int (*lamada)[n],int (*coef)[DEG],int x,int y,int index)
{int temp1[DEG],R[DEG],deone[DEG]={-1};int flag = division(coef[lamada[index][index]],coef[lamada[x][0]],temp1,R);if(!flag) return;remove(lamada,coef,x+1,index+1,temp1,1);remove(lamada,coef,index+1,x+1,deone,1);rl_decline(lamada,coef,y,1,index);
}
//判断左上角元素能否整除所有元素 
//index表示所要处理的第几个左上角元素 以0开始 
int all_division(int (*lamada)[n],int (*coef)[DEG],int index)
{int i,j,flag,prem[DEG],Q[DEG],R[DEG];modify(prem,coef[lamada[index][index]]);for(i=index;i<n;i++){for(j=index;j<n;j++){flag = division(prem,coef[lamada[i][j]],Q,R);if(!flag) break;}}if(i==j) return 1;else return 0;
}
//初等变换使得左上角元素能整除所有元素，本质就是不断降低次数 
void ele_oper(int (*lamada)[n],int (*coef)[DEG],int index)
{int i,j,r_index=index,l_index=index,min = strlenNum(coef[lamada[index][index]]);int prem[DEG],Q[DEG],R[DEG];modify(prem,coef[lamada[index][index]]);for(i=index;i<n;i++){for(j=index;j<n;j++){if(strlenNum(coef[lamada[i][j]]) < min){min = strlenNum(coef[lamada[i][j]]);r_index=i;l_index=j;}}}change(lamada,index+1,r_index+1,1);change(lamada,index+1,l_index+1,0);for(i=index;i<n;i++){if(!division(prem,coef[lamada[i][index]],Q,R) || !division(prem,coef[lamada[index][i]],Q,R)){break;}}if(i<n && !division(prem,coef[lamada[index][i]],Q,R)){rl_decline(lamada,coef,i,1,index);}else if(i<n && !division(prem,coef[lamada[i][index]],Q,R)){rl_decline(lamada,coef,i,0,index);}else if(i==n){for(i=index+1;i<n;i++){for(j=index+1;j<n;j++){if(!division(prem,coef[lamada[i][j]],Q,R))r_index = i;l_index = j;}}all_decline(lamada,coef,r_index,l_index,index);}
}
//将左上角元素变成能整除所有元素 
void middle_temp(int (*lamada)[n],int (*coef)[DEG],int index)
{int len = strlenNum(coef[lamada[index][index]]);int flag = all_division(lamada,coef,index);if(len==1 || flag) return;else {ele_oper(lamada,coef,index);middle_temp(lamada,coef,index);}
}
//开始真正的标准形初等变换
void rel_eleoper(int (*lamada)[n],int (*coef)[DEG],int index)
{int strlen,flag1,flag2,Q[DEG],R[DEG],prem[DEG];modify(prem,coef[lamada[index][index]]);int zeros[DEG] = {0},strp = strlenNum(prem);for(int i=index+1;i<n;i++){strlen = strlenNum(coef[lamada[i][index]]);flag1 = strlen - strp + 1;flag2 = prem[strp-1];int temp[DEG] = {flag1*flag2};division(prem,coef[lamada[i][index]],Q,R);multiple(lamada,coef,temp,i+1,1);remove(lamada,coef,i+1,index+1,Q,1);}for(int i=index+1;i<n;i++){modify(coef[lamada[index][i]],zeros);}
}
//将首项系数化为一结果保存在result中 
void transone(int (*lamada)[n],int (*coef)[DEG],double (*result)[DEG])
{int temp1,temp2;for(int i=0;i<n*n;i++){for(int j=0;j<DEG;j++){result[i][j] = 0;}}for(int i=0;i<n;i++){for(int j=0;j<DEG;j++){temp1 = strlenNum(coef[lamada[i][i]]);temp2 = coef[lamada[i][i]][temp1-1];result[lamada[i][i]][j] = (double)coef[lamada[i][i]][j] / temp2;}}
}int main()
{//系数按从左到右再从上到下排列 int coef[n*n][DEG] = {{1,-1},{-1,2},{0,1},{0,1},{0,0,1},{0,-1},{1,0,1},{-1,1,0,1},{0,0,-1}};//映射lamada矩阵 int lamada[n][n] = {{0,1,2},{3,4,5},{6,7,8}};double result[n*n][DEG];for(int i=0;i<n-1;i++){middle_temp(lamada,coef,i);rel_eleoper(lamada,coef,i);}transone(lamada,coef,result);for(int i=0;i<n;i++){for(int j=0;j<n;j++){for(int k=0;k<DEG;k++){printf("%8.2lf",result[lamada[i][j]][k]);}printf("\n");}}
}