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

void mul_aux(mat_ZZ& X, const mat_ZZ& A, const mat_ZZ& 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 acc, tmp;

    for (i = 1; i <= n; i++) {
        for (j = 1; j <= m; j++) {
            clear(acc);
            for(k = 1; k <= l; k++) {
                mul(tmp, A(i,k), B(k,j));
                add(acc, acc, tmp);
            }
            X(i,j) = acc;
        }
    }
}

void sub(mat_ZZ& X, const mat_ZZ& A, const mat_ZZ& 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));
}

vec_ZZ  GS_reduce(mat_ZZ mat, vec_ZZ vec, ZZ weight)
{
    vec_ZZ          target;
    ZZ_mat<mpz_t>   mat_fp;
    mat_ZZ          mat_exp;
    int             row, col;


    row     =   mat.NumRows();
    col     =   mat.NumCols();
    mat_exp.SetDims(row+1, col+1);

    for(int i=0;i<row;i++)
    {
        for(int j=0;j<col;j++)
        {
            mat_exp[i][j] =   mat[i][j];
        }
    }
    for(int i=0;i<col;i++)
    {
        mat_exp[row][i] =   vec[i];
    }
    mat_exp[row][col]   =   weight;

    mat_fp  =   convert<ZZ_mat<mpz_t> ,mat_ZZ >(mat_exp);
    lllReduction(mat_fp, 0.99,0.51, LM_WRAPPER,FT_DEFAULT, 0, LLL_DEFAULT);
    mat_exp =   convert<mat_ZZ ,ZZ_mat<mpz_t> >(mat_fp);
    target.SetLength(col);
    for (int i=0;i<col;i++)
        target[i]   =   mat_exp[row][i];
//    cout<<target<<endl;
    return target;
}

static
long image(ZZ& det, mat_ZZ& B, mat_ZZ* U, long verbose)
{
   long m = B.NumRows();
   long n = B.NumCols();

   long force_reduce = 1;

   vec_long P;
   P.SetLength(m);

   vec_ZZ D;
   D.SetLength(m+1);
   D[0] = 1;

   vec_vec_ZZ lam;

   lam.SetLength(m);

   long j;
   for (j = 1; j <= m; j++)
      lam(j).SetLength(m);

   if (U) ident(*U, m);

   long s = 0;

   long k = 1;
   long max_k = 0;


   while (k <= m) {
      if (k > max_k) {
         IncrementalGS(B, P, D, lam, s, k);
         max_k = k;
      }

      if (k == 1) {
         force_reduce = 1;
         k++;
         continue;
      }

      if (force_reduce)
         for (j = k-1; j >= 1; j--)
            reduce(k, j, B, P, D, lam, U);

      if (P(k-1) != 0 && P(k) == 0) {
         force_reduce = swap(k, B, P, D, lam, U, max_k, verbose);
         k--;
      }
      else {
         force_reduce = 1;
         k++;
      }
   }

   det = D[s];
   return s;
}

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

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

    mat_ZZ 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& X, const mat_ZZ& 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 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);
    }
}

NTL_START_IMPL


void add(mat_ZZ& X, const mat_ZZ& A, const mat_ZZ& 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));  
}  

void conv(mat_ZZ_p& x, const mat_ZZ& a)
{
    long n = a.NumRows();
    long m = a.NumCols();
    long i;

    x.SetDims(n, m);
    for (i = 0; i < n; i++)
        conv(x[i], a[i]);
}

static
long MaxBits(const mat_ZZ& A)
{
    long m = 0;
    long i, j;
    for (i = 0; i < A.NumRows(); i++)
        for (j = 0; j < A.NumCols(); j++)
            m = max(m, NumBits(A[i][j]));

    return m;
}

void mul(mat_ZZ& X, const mat_ZZ& A, long b)
{
    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);
}

static
void ExactDiv(mat_ZZ& x, const ZZ& d)
{
    long n = x.NumRows();
    long m = x.NumCols();

    long i, j;

    for (i = 0; i < n; i++)
        for (j = 0; j < m; j++)
            if (!divide(x[i][j], x[i][j], d))
                Error("inexact division");
}

void negate(mat_ZZ& X, const mat_ZZ& 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));
}

static
void ComputeGS(mat_ZZ& B, xdouble **B1, xdouble **mu, xdouble *b, 
               xdouble *c, long k, xdouble bound, long st, xdouble *buf)
{
   long n = B.NumCols();
   long i, j;
   xdouble s, t1, y, t;
   ZZ T1;

   xdouble *mu_k = mu[k];

   if (st < k) {
      for (i = 1; i < st; i++)
         buf[i] = mu_k[i]*c[i];
   }

   for (j = st; j <= k-1; j++) {
      if (b[k]*b[j] < NTL_FDOUBLE_PRECISION*NTL_FDOUBLE_PRECISION) {
         double z = 0;
         xdouble *B1_k = B1[k];
         xdouble *B1_j = B1[j];

         for (i = 1; i <= n; i++)
            z += B1_k[i].x * B1_j[i].x;

         s = z;
      }
      else {
         s = InnerProduct(B1[k], B1[j], n);
   
         if (s*s <= b[k]*b[j]/bound) {
            InnerProduct(T1, B(k), B(j));
            conv(s, T1);
         }
      }

      xdouble *mu_j = mu[j];

      t1 = 0;
      for (i = 1; i <= j-1; i++)
         MulAdd(t1, t1, mu_j[i], buf[i]);

      mu_k[j] = (buf[j] = (s - t1))/c[j];
   }

   s = 0;
   for (j = 1; j <= k-1; j++)
      MulAdd(s, s, mu_k[j], buf[j]);

   c[k] = b[k] - s;
}

static
void IncrementalGS(mat_ZZ& B, vec_long& P, vec_ZZ& D, vec_vec_ZZ& lam,
                   long& s, long k)
{
   long n = B.NumCols();
   long m = B.NumRows();

   static ZZ u, t1, t2;

   long i, j;

   for (j = 1; j <= k-1; j++) {
      long posj = P(j);
      if (posj == 0) continue;

      InnerProduct(u, B(k), B(j));
      for (i = 1; i <= posj-1; i++) {
         mul(t1, D[i], u);
         mul(t2, lam(k)(i), lam(j)(i));
         sub(t1, t1, t2);
         div(t1, t1, D[i-1]);
         u = t1;
      }

      lam(k)(posj) = u;
   }

   InnerProduct(u, B(k), B(k));
   for (i = 1; i <= s; i++) {
      mul(t1, D[i], u);
      mul(t2, lam(k)(i), lam(k)(i));
      sub(t1, t1, t2);
      div(t1, t1, D[i-1]);
      u = t1;
   }

   if (u == 0) {
      P(k) = 0;
   }
   else {
      s++;
      P(k) = s;
      D[s] = u;
   }
}

long IsIdent(const mat_ZZ& A, long n)
{
    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 (!IsOne(A(i, j))) return 0;
            }

    return 1;
}