NTL密码算法开源库——模二整数上的矩阵（mat_GF2）

2021SC@SDUSC

矩阵运算

具体代码

#include <NTL/matrix.h>
#include <NTL/vec_vec_GF2.h>


typedef Mat<GF2> mat_GF2; // backward compatibility


void conv(mat_GF2& X, const vec_vec_GF2& A);
mat_GF2 to_mat_GF2(const vec_vec_GF2& A);
// convert a vector of vec_GF2's to a matrix


// procedural arithmetic routines:

void add(mat_GF2& X, const mat_GF2& A, const mat_GF2& B);
// X = A + B

void sub(mat_GF2& X, const mat_GF2& A, const mat_GF2& B);
// X = A - B = A + B

void negate(mat_GF2& X, const mat_GF2& A);
// X = -A = A 

void mul(mat_GF2& X, const mat_GF2& A, const mat_GF2& B);
// X = A * B

void mul(vec_GF2& x, const mat_GF2& A, const vec_GF2& b);
// x = A * b

void mul(vec_GF2& x, const vec_GF2& a, const mat_GF2& B);
// x = a * B


void mul(mat_GF2& X, const mat_GF2& A, GF2 b);
void mul(mat_GF2& X, const mat_GF2& A, long b);
// X = A * b

void mul(mat_GF2& X, GF2 a, const mat_GF2& B);
void mul(mat_GF2& X, long a, const mat_GF2& B);
// X = a * B

void add(mat_GF2& X, const mat_GF2& A, const mat_GF2& B)  
{  
   long n = A.NumRows();  
   long m = A.NumCols();  
  
   if (B.NumRows() != n || B.NumCols() != m)   
      LogicError("matrix add: dimension mismatch");  
  
   X.SetDims(n, m);  

   long mw = (m + NTL_BITS_PER_LONG - 1)/NTL_BITS_PER_LONG;
  
   long i;  
   for (i = 0; i < n; i++) {
      _ntl_ulong *xp = X[i].rep.elts();
      const _ntl_ulong *ap = A[i].rep.elts();
      const _ntl_ulong *bp = B[i].rep.elts();
      long j;
      for (j = 0; j < mw; j++)
         xp[j] = ap[j] ^ bp[j];
   }
}  
void mul(vec_GF2& x, const mat_GF2& A, const vec_GF2& b)  
{  
   if (&b == &x || A.alias(x)) {
      vec_GF2 tmp;
      mul_aux(tmp, A, b);
      x = tmp;
   }
   else
      mul_aux(x, A, b);
}  
void mul(vec_GF2& x, const vec_GF2& a, const mat_GF2& B)
{
   if (&a == &x || B.alias(x)) {
      vec_GF2 tmp;
      mul_aux(tmp, a, B);
      x = tmp;
   }
   else
      mul_aux(x, a, B);
}
void ident(mat_GF2& X, long n)  
{  
   X.SetDims(n, n);  
   clear(X);
   long i;  
  
   for (i = 0; i < n; i++)  
      X.put(i, i, to_GF2(1));
} 


void determinant(ref_GF2 d, const mat_GF2& M_in)
{
   long k, n;
   long i, j;
   long pos;

   n = M_in.NumRows();

   if (M_in.NumCols() != n)
      LogicError("determinant: nonsquare matrix");

   if (n == 0) {
      set(d);
      return;
   }

   mat_GF2 M;

   M = M_in;

   long wn = (n + NTL_BITS_PER_LONG - 1)/NTL_BITS_PER_LONG;

   for (k = 0; k < n; k++) {
      long wk = k/NTL_BITS_PER_LONG;
      long bk = k - wk*NTL_BITS_PER_LONG;
      _ntl_ulong k_mask = 1UL << bk;

      pos = -1;
      for (i = k; i < n; i++) {
         if (M[i].rep.elts()[wk] & k_mask) {
            pos = i;
            break;
         }
      }

      if (pos != -1) {
         if (k != pos) {
            swap(M[pos], M[k]);
         }


         _ntl_ulong *y = M[k].rep.elts();

         for (i = k+1; i < n; i++) {
            // M[i] = M[i] + M[k]*M[i,k]

            if (M[i].rep.elts()[wk] & k_mask) {
               _ntl_ulong *x = M[i].rep.elts();

               for (j = wk; j < wn; j++)
                  x[j] ^= y[j];
            }

         }
      }
      else {
         clear(d);
         return;
      }
   }

   set(d);
   return;
}

具体数学原理见@回首，阑珊

高斯消元

1.简介
高斯消元是一个用于求方程组的解的算法。在线性代数中非常重要。一般而言，其复杂度为O ( n 3 )可以承受200 200200及以下的数据范围(当然有的题目时限是10000 m s 10000ms10000ms什么的，特殊情况特殊对待)。
高斯消元法(Gaussian elimination)是求解线性方阵组的一种算法，它也可用来求矩阵的秩，以及求可逆方阵的逆矩阵。它通过逐步消除未知数来将原始线性系统转化为另一个更简单的等价的系统。它的实质是通过初等行变化(Elementary row operations)，将线性方程组的增广矩阵转化为行阶梯矩阵(row echelon form).
2.算法理解
考虑一个二元一次方程组怎么解。

应用：整数线性方程组
具体代码

int GCD(int x,int y){return x%y==0?y:GCD(y,x%y);}
int LCM(int x,int y){return x/GCD(x,y)*y;}
int a[MAXN][MAXN],x[MAXN];
bool frex[MAXN];
int Gauss(int equ,int var){
	for(int i=0;i<=var;i++) x[i]=0,frex[i]=1;
	//step1
	int col=0,row;
	for(row=0;row<equ&&col<var;row++,col++){
		//step1.1
		int mx_r=row;
		for(int i=row+1;i<equ;i++)
			if(abs(a[i][col])>abs(a[mx_r][col]))
				mx_r=i;
		if(mx_r!=row)
			for(int j=col;j<=var;j++)
				swap(a[row][j],a[mx_r][j]);
		if(a[row][col]==0){row--;continue;}
		//step1.2
		for(int i=row+1;i<equ;i++)
			if(a[i][col]!=0){
				int lcm=LCM(abs(a[i][col]),abs(a[row][col]));
				int ta=lcm/abs(a[i][col]),tb=lcm/abs(a[row][col]);//为了避免精度损失，要用lcm
				if(a[i][col]*a[row][col]<0) tb=-tb;
				for(int j=col;j<=var;j++)
					a[i][j]=a[i][j]*ta-a[row][j]*tb;
			}
	}
	//step2
	for(int i=row;i<equ;i++)
		if(a[i][col]!=0)
			return -1;//step2.1
	if(var>row) return var-row;//step2.2
	//step2.3
	for(int i=var-1;i>=0;i--){
		int tmp=a[i][var];
		for(int j=i+1;j<var;j++)
			tmp-=x[j]*a[i][j];
		if(tmp%a[i][i]!=0) return -2;//如果无法整除，则无整数解
		x[i]=tmp/a[i][i];
	}return 0;
}