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

/**
 * 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);
}

Mat<ZZ> CreateLPerp(Mat<ZZ> const &A, int n, int m, long q, double sigma) {
	Mat<ZZ> LPerp;
	LPerp.SetDims(2 * n + m, n + m);

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

	ZZ factor = (2 * conv<ZZ>(sigma));

	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n + m; j++) {
			if (j < m)
				LPerp[i][j] = -A[i][j] % q;
			if (j - i == m)
				LPerp[i][j] = factor;
		}
	}

	for (int i = n; i < 2 * n + m; i++) {
		for (int j = 0; j < n + m; j++)
			if (j < m)
				LPerp[i][j] = qI[i - n][j];
			else
				LPerp[i][j] = factor * qI[i - n][j];
	}

	return LPerp;
}

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 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 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> PadMatrix(Mat<ZZ> A, long q, int m, int n) {
	Mat<ZZ> res;
	Mat<ZZ> qI;
	res.SetDims(m + n, m);

	A = transpose(A);

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

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

	for (int i = n; i < m + n; i++) {
		for (int j = 0; j < m; j++) {
			res[i][j] = qI[i - n][j];
		}
	}
	return res;
}

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;
}

int main()
{
	int iA[MS][MS] = {	{5,2,1},
						{3,0,-4},
						{2,1,6}};

	Mat<ZZ> A;
	long i, j;
	ZZ d;

	A.SetDims( MS, MS );

	for (i = 0; i < MS; i++)
		for (j = 0; j < MS; j++)
			A[i][j] = iA[i][j];

	cout << A << endl;

	long deterministic=0;
	d = determinant(A, deterministic);

	cout << d << endl;

}

/**
 * GenerateSamples
 * \brief generates the requested number of samples and returns them in matrix
 *        form
 *
 * @param n length of each sample vector
 * @param q modulus
 * @param m number of samples
 * @param sigma for discrete gaussian
 * @param tailfactor for discrete gaussian
 * @param number_of_rects used in sampling algorithm (trade off speed/memory)
 * @param precision used in sampling algorithm
 */
tuple<Mat<ZZ>, Vec<ZZ>, Vec<ZZ>>
GenerateSamples(int n, long q, int m, double sigma, int tailfactor,
				int number_of_rects, int precision, int flag_binary,
				int flag_trinary, int flag_binary_secret, int flag_binary_lwe,
				int flag_binary_sis) {
	Mat<ZZ> A;
	Vec<ZZ> t;
	Vec<ZZ> s;
	Vec<ZZ> e;
	Vec<ZZ> As; // solution vector -EK

	Vec<ZZ> u; // to test BinaryA
	Vec<ZZ> C1; // to test BinaryA
	Mat<ZZ> Abin;

	u.SetLength(m);
	C1.SetLength(n);
	Abin.SetDims(n, m);

	A.SetDims(m, n); // m rows x n cols
	t.SetLength(m);
	As.SetLength(m);
	// generate LWE secret s
	if (flag_binary_secret) {
		s = RandomVec(n, 2);
	} else {
		s = RandomVec(n, q);
	}

	// generate error vector e
	if (flag_binary) {
		random_device rd;
		uniform_int_distribution<int> dist(0, 1);
		e.SetLength(m);
		for (ZZ *it = e.begin(); it != e.end(); ++it) {
			*it = ZZ(dist(rd));
		}
		// a uniform, trinary error (from {-1, 0, 1}) is used
		// in "Efficient Identity-Based Encryption over NTRU Lattices"
		// by Ducas, Lyubashevsky and Prest
	} else if (flag_trinary) {
		random_device rd;
		uniform_int_distribution<int> dist(-1, 1);
		e.SetLength(m);
		for (ZZ *it = e.begin(); it != e.end(); ++it) {
			*it = ZZ(dist(rd));
		}
	} else {
		Ziggurat sampler(number_of_rects, sigma, tailfactor, precision);
		e.SetLength(m);
		for (ZZ *it = e.begin(); it != e.end(); ++it) {
			*it = sampler.sample();
		}
	}

	// generate LWE samples
	if (flag_binary_sis) {
		//-------For Binary-SIS---------
		u = RandomVec(m, 2);

		for (int i = 0; i < n; i++) {
			Abin[i] = RandomVec(m, 2);
			InnerProduct(C1[i], Abin[i], u);
			C1[i] = C1[i] % q;
		}
		return make_tuple(A, t, As);
	} else if (flag_binary_lwe) {
		//-------For Binary-LWE---------
		for (int i = 0; i < m; ++i) {
			A[i] = RandomVec(n, 2);
			InnerProduct(t[i], A[i], s);
			As[i] = t[i] % q;
			t[i] += e[i];
			t[i] = t[i] % q;
		}
		return make_tuple(A, t, As);
	} else {
		for (int i = 0; i < m; ++i) {
			A[i] = RandomVec(n, q);
			InnerProduct(t[i], A[i], s);
			As[i] = t[i] % q;
			t[i] += e[i];
			t[i] = t[i] % q;
		}
		return make_tuple(A, t, As);
	}
}

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;
//.........这里部分代码省略.........