C++ mat_zz_pE类代码示例

void power(mat_zz_pE& X, const mat_zz_pE& A, const ZZ& e)
{
   if (A.NumRows() != A.NumCols()) Error("power: non-square matrix");

   if (e == 0) {
      ident(X, A.NumRows());
      return;
   }

   mat_zz_pE T1, T2;
   long i, k;

   k = NumBits(e);
   T1 = A;

   for (i = k-2; i >= 0; i--) {
      sqr(T2, T1);
      if (bit(e, i))
         mul(T1, T2, A);
      else
         T1 = T2;
   }

   if (e < 0)
      inv(X, T1);
   else
      X = T1;
}

void transpose(mat_zz_pE& X, const mat_zz_pE& A)
{
   long n = A.NumRows();
   long m = A.NumCols();

   long i, j;

   if (&X == & A) {
      if (n == m)
         for (i = 1; i <= n; i++)
            for (j = i+1; j <= n; j++)
               swap(X(i, j), X(j, i));
      else {
         mat_zz_pE tmp;
         tmp.SetDims(m, n);
         for (i = 1; i <= n; i++)
            for (j = 1; j <= m; j++)
               tmp(j, i) = A(i, j);
         X.kill();
         X = tmp;
      }
   }
   else {
      X.SetDims(m, n);
      for (i = 1; i <= n; i++)
         for (j = 1; j <= m; j++)
            X(j, i) = A(i, j);
   }
}

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())  
      Error("matrix mul: dimension mismatch");  
  
   X.SetDims(n, m);  
  
   long i, j, k;  
   zz_pX acc, tmp;  
  
   for (i = 1; i <= n; i++) {  
      for (j = 1; j <= m; j++) {  
         clear(acc);  
         for(k = 1; k <= l; k++) {  
            mul(tmp, rep(A(i,k)), rep(B(k,j)));  
            add(acc, acc, tmp);  
         }  
         conv(X(i,j), acc);  
      }  
   }  
}  

void negate(mat_zz_pE& X, const mat_zz_pE& A)  
{  
   long n = A.NumRows();  
   long m = A.NumCols();  
  
  
   X.SetDims(n, m);  
  
   long i, j;  
   for (i = 1; i <= n; i++)  
      for (j = 1; j <= m; j++)  
         negate(X(i,j), A(i,j));  
}  

void mul(mat_zz_pE& X, const mat_zz_pE& A, const zz_pE& b_in)
{
   zz_pE b = b_in;
   long n = A.NumRows();
   long m = A.NumCols();

   X.SetDims(n, m);

   long i, j;
   for (i = 0; i < n; i++)
      for (j = 0; j < m; j++)
         mul(X[i][j], A[i][j], b);
}

void sub(mat_zz_pE& X, const mat_zz_pE& A, const mat_zz_pE& B)  
{  
   long n = A.NumRows();  
   long m = A.NumCols();  
  
   if (B.NumRows() != n || B.NumCols() != m)  
      Error("matrix sub: dimension mismatch");  
  
   X.SetDims(n, m);  
  
   long i, j;  
   for (i = 1; i <= n; i++)  
      for (j = 1; j <= m; j++)  
         sub(X(i,j), A(i,j), B(i,j));  
}  

void convert(vector< vector<ZZX> >& X, const mat_zz_pE& A)
{
   long n = A.NumRows();
   X.resize(n);
   for (long i = 0; i < n; i++)
      convert(X[i], A[i]);
}

void clear(mat_zz_pE& x)
{
   long n = x.NumRows();
   long i;
   for (i = 0; i < n; i++)
      clear(x[i]);
}

long IsDiag(const mat_zz_pE& A, long n, const zz_pE& d)
{
   if (A.NumRows() != n || A.NumCols() != n)
      return 0;

   long i, j;

   for (i = 1; i <= n; i++)
      for (j = 1; j <= n; j++)
         if (i != j) {
            if (!IsZero(A(i, j))) return 0;
         }
         else {
            if (A(i, j) != d) return 0;
         }

   return 1;
}

NTL_START_IMPL

  
void add(mat_zz_pE& X, const mat_zz_pE& A, const mat_zz_pE& 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 i, j;  
   for (i = 1; i <= n; i++)   
      for (j = 1; j <= m; j++)  
         add(X(i,j), A(i,j), B(i,j));  
}  

long IsZero(const mat_zz_pE& a)
{
   long n = a.NumRows();
   long i;

   for (i = 0; i < n; i++)
      if (!IsZero(a[i]))
         return 0;

   return 1;
}

// prime power solver
// A is an n x n matrix, we compute its inverse mod p^r. An error is raised
// if A is not inverible mod p. zz_p::modulus() is assumed to be p^r, for
// p prime, r >= 1. Also zz_pE::modulus() is assumed to be initialized.
void ppInvert(mat_zz_pE& X, const mat_zz_pE& A, long p, long r)
{
  if (r == 1) { // use native inversion from NTL
    inv(X, A);    // X = A^{-1}
    return;
  }

  // begin by inverting A modulo p

  // convert to ZZX for a safe transaltion to mod-p objects
  vector< vector<ZZX> > tmp;
  convert(tmp, A);
  { // open a new block for mod-p computation
  ZZX G;
  convert(G, zz_pE::modulus());
  zz_pBak bak_pr; bak_pr.save(); // backup the mod-p^r moduli
  zz_pEBak bak_prE; bak_prE.save();
  zz_p::init(p);   // Set the mod-p moduli
  zz_pE::init(conv<zz_pX>(G));

  mat_zz_pE A1, Inv1;
  convert(A1, tmp);   // Recover A as a mat_zz_pE object modulo p
  inv(Inv1, A1);      // Inv1 = A^{-1} (mod p)
  convert(tmp, Inv1); // convert to ZZX for transaltion to a mod-p^r object
  } // mod-p^r moduli restored on desctuction of bak_pr and bak_prE
  mat_zz_pE XX;
  convert(XX, tmp); // XX = A^{-1} (mod p)

  // Now lift the solution modulo p^r

  // Compute the "correction factor" Z, s.t. XX*A = I - p*Z (mod p^r)
  long n = A.NumRows();
  const mat_zz_pE I = ident_mat_zz_pE(n); // identity matrix
  mat_zz_pE Z = I - XX*A;

  convert(tmp, Z);  // Conver to ZZX to divide by p
  for (long i=0; i<n; i++) for (long j=0; j<n; j++) tmp[i][j] /= p;
  convert(Z, tmp);  // convert back to a mod-p^r object

  // The inverse of A is ( I+(pZ)+(pZ)^2+...+(pZ)^{r-1} )*XX (mod p^r). We use
  // O(log r) products to copmute it as (I+pZ)* (I+(pZ)^2)* (I+(pZ)^4)*...* XX

  long e = NextPowerOfTwo(r); // 2^e is smallest power of two >= r

  Z *= p;                 // = pZ
  mat_zz_pE prod = I + Z; // = I + pZ
  for (long i=1; i<e; i++) {
    sqr(Z, Z);     // = (pZ)^{2^i}
    prod *= (I+Z); // = sum_{j=0}^{2^{i+1}-1} (pZ)^j
  }
  mul(X, prod, XX); // X = A^{-1} mod p^r
  assert(X*A == I);
}

void ident(mat_zz_pE& X, long n)  
{  
   X.SetDims(n, n);  
   long i, j;  
  
   for (i = 1; i <= n; i++)  
      for (j = 1; j <= n; j++)  
         if (i == j)  
            set(X(i, j));  
         else  
            clear(X(i, j));  
} 

void diag(mat_zz_pE& X, long n, const zz_pE& d_in)  
{  
   zz_pE d = d_in;
   X.SetDims(n, n);  
   long i, j;  
  
   for (i = 1; i <= n; i++)  
      for (j = 1; j <= n; j++)  
         if (i == j)  
            X(i, j) = d;  
         else  
            clear(X(i, j));  
} 

void buildLinPolyMatrix(mat_zz_pE& M, long p)
{
   long d = zz_pE::degree();

   M.SetDims(d, d);

   for (long j = 0; j < d; j++) 
      conv(M[0][j], zz_pX(j, 1));

   for (long i = 1; i < d; i++)
      for (long j = 0; j < d; j++)
         M[i][j] = power(M[i-1][j], p);
}