C++ Mat::NumCols方法代码示例

/**
 * muGSO
 * \brief on input B returnes a lower-triangular matrix mu such that B = mu*Q,
 *        with <b_i^*, b_i^*> / ||b_i^*|| = 1 on the main diagonal,
 *        Q - matrix with orthonormal rows.
 *
 * @params B the basis
 */
matrix<double> muGSO(matrix<long> const &B_std) {
	Mat<ZZ> B = to_ntl<ZZ>(B_std);
	Mat<RR> B_star;
	Mat<RR> muGSO;
	muGSO.SetDims(B.NumRows(), B.NumCols());
	B_star.SetDims(B.NumRows(), B.NumCols());
	if (B.NumRows() > 0) {
		for (long i = 0; i < B.NumRows(); ++i) {
			B_star[i] = conv<Vec<RR>>(B[i]);
			for (long j = 0; j < i; ++j) {
				RR num; num = 0;
				RR denom; denom = 0;

				InnerProduct(num, B_star[j], conv<Vec<RR>>(B[i]));
				InnerProduct(denom, B_star[j], B_star[j]);

				muGSO[i][j] = num / denom;
				B_star[i] = B_star[i] - muGSO[i][j] * B_star[j];
			}
			// InnerProduct(muGSO[i][i], B_star[i], B_star[i]);
			muGSO[i][i] = 1;
		}
	}
	return to_stl<double>(muGSO);
}

static
void mul_aux(Mat<zz_pE>& X, const Mat<zz_pE>& A, const Mat<zz_pE>& B)  
{  
   long n = A.NumRows();  
   long l = A.NumCols();  
   long m = B.NumCols();  
  
   if (l != B.NumRows())  
      LogicError("matrix mul: dimension mismatch");  
  
   X.SetDims(n, m);  


   zz_pContext zz_p_context;
   zz_p_context.save();
   zz_pEContext zz_pE_context;
   zz_pE_context.save();
   double sz = zz_pE_SizeInWords();

   bool seq = (double(n)*double(l)*double(m)*sz*sz < PAR_THRESH);

   NTL_GEXEC_RANGE(seq, m, first, last)
   NTL_IMPORT(n)
   NTL_IMPORT(l)
   NTL_IMPORT(m)

   zz_p_context.restore();
   zz_pE_context.restore();

   long i, j, k;  
   zz_pX acc, tmp;  

   Vec<zz_pE> B_col;
   B_col.SetLength(l);

   for (j = first; j < last; j++) {
      for (k = 0; k < l; k++) B_col[k] = B[k][j];

      for (i = 0; i < n; i++) {
         clear(acc);
         for (k = 0; k < l; k++) {
            mul(tmp, rep(A[i][k]), rep(B_col[k]));
            add(acc, acc, tmp);
         }
         conv(X[i][j], acc);
      }
   }

   NTL_GEXEC_RANGE_END
}  

static
void RawConvertTranspose(Mat<zz_p>& X, const Mat<MatPrime_residue_t>& A)
{
   long n = A.NumRows();
   long m = A.NumCols();

   X.SetDims(m, n);
   for (long i = 0; i < n; i++) {
      const MatPrime_residue_t *Ai = A[i].elts();
      for (long j = 0; j < m; j++)
         X[j][i] = Ai[j];
   }
} 

static
void RawConvert(Mat<MatPrime_residue_t>& X, const Mat<zz_p>& A)
{
   long n = A.NumRows();
   long m = A.NumCols();

   X.SetDims(n, m);
   for (long i = 0; i < n; i++) {
      const zz_p *Ai = A[i].elts();
      MatPrime_residue_t *Xi = X[i].elts();
      for (long j = 0; j < m; j++)
         Xi[j] = rep(Ai[j]);
   }
} 

Mat<ZZ> Kernel(Mat<ZZ> U, long r, long q) {
	Mat<ZZ> ker;
	long m = U.NumCols();
	// ker.SetDims(m-r, m); //without q-ary part
	ker.SetDims(2 * m - r, m);
	for (int i = 0; i < m - r; i++)
		for (int j = 0; j < m; j++)
			ker[i][j] = U[i][j] % q;

	Mat<ZZ> qI;
	diag(qI, conv<long>(m), conv<ZZ>(q));

	for (int i = m - r; i < 2 * m - r; i++) {
		for (int j = 0; j < m; j++) {
			ker[i][j] = qI[i - m + r][j];
		}
	}

	return ker;
}

/**
 * GSO
 * \brief returns the Gram-Schmidt orthogonalized matrix
 *
 * @params B the basis
 */
matrix<double> GSO(matrix<long> const &B_std) {
	Mat<ZZ> B = to_ntl<ZZ>(B_std);
	Mat<RR> B_star;
	B_star.SetDims(B.NumRows(), B.NumCols());
	if (B.NumRows() > 0) {
		B_star[0] = conv<Vec<RR>>(B[0]);
		for (int i = 1; i < B.NumRows(); ++i) {
			B_star[i] = conv<Vec<RR>>(B[i]);
			for (int j = 0; j < i; ++j) {
				RR num; num = 0;
				RR denom; denom = 0;
				RR division; division = 0;

				InnerProduct(num, B_star[j], conv<Vec<RR>>(B[i]));
				InnerProduct(denom, B_star[j], B_star[j]);

				division = num / denom;

				B_star[i] = B_star[i] - division * B_star[j];
			}
		}
	}
	return to_stl<double>(B_star);
}

void inv(zz_pE& d, Mat<zz_pE>& X, const Mat<zz_pE>& A)
{
   long n = A.NumRows();

   if (A.NumCols() != n)
      LogicError("inv: nonsquare matrix");

   if (n == 0) {
      set(d);
      X.SetDims(0, 0);
      return;
   }

   const zz_pXModulus& G = zz_pE::modulus();

   zz_pX t1, t2;
   zz_pX pivot;
   zz_pX pivot_inv;

   Vec< Vec<zz_pX> > M;
   // scratch space

   M.SetLength(n);
   for (long i = 0; i < n; i++) {
      M[i].SetLength(n);
      for (long j = 0; j < n; j++) {
         M[i][j].SetMaxLength(2*deg(G)-1);
         M[i][j] = rep(A[i][j]);
      }
   }

   zz_pX det;
   det = 1;


   Vec<long> P;
   P.SetLength(n);
   for (long k = 0; k < n; k++) P[k] = k;
   // records swap operations
   

   zz_pContext zz_p_context;
   zz_p_context.save();
   double sz = zz_pE_SizeInWords();

   bool seq = double(n)*double(n)*sz*sz < PAR_THRESH;

   bool pivoting = false;

   for (long k = 0; k < n; k++) {

      long pos = -1;

      for (long i = k; i < n; i++) {
         rem(pivot, M[i][k], G);
         if (pivot != 0) {
            InvMod(pivot_inv, pivot, G);
            pos = i;
            break;
         }
      }

      if (pos != -1) {
         if (k != pos) {
            swap(M[pos], M[k]);
            negate(det, det); 
            P[k] = pos;
            pivoting = true;
         }

         MulMod(det, det, pivot, G);

         {
            // multiply row k by pivot_inv
            zz_pX *y = &M[k][0];
            for (long j = 0; j < n; j++) {
               rem(t2, y[j], G);
               MulMod(y[j], t2, pivot_inv, G);
            }
            y[k] = pivot_inv;
         }


         NTL_GEXEC_RANGE(seq, n, first, last)  
         NTL_IMPORT(n)
         NTL_IMPORT(k)

         zz_p_context.restore();

         zz_pX *y = &M[k][0]; 
         zz_pX t1, t2;

         for (long i = first; i < last; i++) {
            if (i == k) continue; // skip row k

            zz_pX *x = &M[i][0]; 
            rem(t1, x[k], G);
            negate(t1, t1); 
            x[k] = 0;
            if (t1 == 0) continue;
//.........这里部分代码省略.........