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

long CRT(mat_ZZ& gg, ZZ& a, const mat_zz_p& G)
{
    long n = gg.NumRows();
    long m = gg.NumCols();

    if (G.NumRows() != n || G.NumCols() != m)
        Error("CRT: dimension mismatch");

    long p = zz_p::modulus();

    ZZ new_a;
    mul(new_a, a, p);

    long a_inv;
    a_inv = rem(a, p);
    a_inv = InvMod(a_inv, p);

    long p1;
    p1 = p >> 1;

    ZZ a1;
    RightShift(a1, a, 1);

    long p_odd = (p & 1);

    long modified = 0;

    long h;

    ZZ g;
    long i, j;

    for (i = 0; i < n; i++) {
        for (j = 0; j < m; j++) {
            if (!CRTInRange(gg[i][j], a)) {
                modified = 1;
                rem(g, gg[i][j], a);
                if (g > a1) sub(g, g, a);
            }
            else
                g = gg[i][j];

            h = rem(g, p);
            h = SubMod(rep(G[i][j]), h, p);
            h = MulMod(h, a_inv, p);
            if (h > p1)
                h = h - p;

            if (h != 0) {
                modified = 1;

                if (!p_odd && g > 0 && (h == p1))
                    MulSubFrom(g, a, h);
                else
                    MulAddTo(g, a, h);

            }

            gg[i][j] = g;
        }
    }

    a = new_a;

    return modified;

}

void solve(ZZ& d_out, vec_ZZ& x_out,
           const mat_ZZ& A, const vec_ZZ& b,
           long deterministic)
{
    long n = A.NumRows();

    if (A.NumCols() != n)
        Error("solve: nonsquare matrix");

    if (b.length() != n)
        Error("solve: dimension mismatch");

    if (n == 0) {
        set(d_out);
        x_out.SetLength(0);
        return;
    }

    zz_pBak zbak;
    zbak.save();

    ZZ_pBak Zbak;
    Zbak.save();

    vec_ZZ x(INIT_SIZE, n);
    ZZ d, d1;

    ZZ d_prod, x_prod;
    set(d_prod);
    set(x_prod);

    long d_instable = 1;
    long x_instable = 1;

    long check = 0;

    long gp_cnt = 0;

    vec_ZZ y, b1;

    long i;
    long bound = 2+DetBound(A);

    for (i = 0; ; i++) {
        if ((check || IsZero(d)) && !d_instable) {
            if (NumBits(d_prod) > bound) {
                break;
            }
            else if (!deterministic &&
                     bound > 1000 && NumBits(d_prod) < 0.25*bound) {

                ZZ P;

                long plen = 90 + NumBits(max(bound, NumBits(d)));
                GenPrime(P, plen, 90 + 2*NumBits(gp_cnt++));

                ZZ_p::init(P);

                mat_ZZ_p AA;
                conv(AA, A);

                ZZ_p dd;
                determinant(dd, AA);

                if (CRT(d, d_prod, rep(dd), P))
                    d_instable = 1;
                else
                    break;
            }
        }


        zz_p::FFTInit(i);
        long p = zz_p::modulus();

        mat_zz_p AA;
        conv(AA, A);

        if (!check) {
            vec_zz_p bb, xx;
            conv(bb, b);

            zz_p dd;

            solve(dd, xx, AA, bb);

            d_instable = CRT(d, d_prod, rep(dd), p);
            if (!IsZero(dd)) {
                mul(xx, xx, dd);
                x_instable = CRT(x, x_prod, xx);
            }
            else
                x_instable = 1;

            if (!d_instable && !x_instable) {
                mul(y, x, A);
                mul(b1, b, d);
                if (y == b1) {
                    d1 = d;
                    check = 1;
//.........这里部分代码省略.........

void solve1(ZZ& d_out, vec_ZZ& x_out, const mat_ZZ& A, const vec_ZZ& b)
{
    long n = A.NumRows();

    if (A.NumCols() != n)
        Error("solve1: nonsquare matrix");

    if (b.length() != n)
        Error("solve1: dimension mismatch");

    if (n == 0) {
        set(d_out);
        x_out.SetLength(0);
        return;
    }

    ZZ num_bound, den_bound;

    hadamard(num_bound, den_bound, A, b);

    if (den_bound == 0) {
        clear(d_out);
        return;
    }

    zz_pBak zbak;
    zbak.save();

    long i;
    long j;

    ZZ prod;
    prod = 1;

    mat_zz_p B;


    for (i = 0; ; i++) {
        zz_p::FFTInit(i);

        mat_zz_p AA, BB;
        zz_p dd;

        conv(AA, A);
        inv(dd, BB, AA);

        if (dd != 0) {
            transpose(B, BB);
            break;
        }

        mul(prod, prod, zz_p::modulus());

        if (prod > den_bound) {
            d_out = 0;
            return;
        }
    }

    long max_A_len = MaxBits(A);

    long use_double_mul1 = 0;
    long use_double_mul2 = 0;
    long double_limit = 0;

    if (max_A_len + NTL_SP_NBITS + NumBits(n) <= NTL_DOUBLE_PRECISION-1)
        use_double_mul1 = 1;

    if (!use_double_mul1 && max_A_len+NTL_SP_NBITS+2 <= NTL_DOUBLE_PRECISION-1) {
        use_double_mul2 = 1;
        double_limit = (1L << (NTL_DOUBLE_PRECISION-1-max_A_len-NTL_SP_NBITS));
    }

    long use_long_mul1 = 0;
    long use_long_mul2 = 0;
    long long_limit = 0;

    if (max_A_len + NTL_SP_NBITS + NumBits(n) <= NTL_BITS_PER_LONG-1)
        use_long_mul1 = 1;

    if (!use_long_mul1 && max_A_len+NTL_SP_NBITS+2 <= NTL_BITS_PER_LONG-1) {
        use_long_mul2 = 1;
        long_limit = (1L << (NTL_BITS_PER_LONG-1-max_A_len-NTL_SP_NBITS));
    }



    if (use_double_mul1 && use_long_mul1)
        use_long_mul1 = 0;
    else if (use_double_mul1 && use_long_mul2)
        use_long_mul2 = 0;
    else if (use_double_mul2 && use_long_mul1)
        use_double_mul2 = 0;
    else if (use_double_mul2 && use_long_mul2) {
        if (long_limit > double_limit)
            use_double_mul2 = 0;
        else
            use_long_mul2 = 0;
    }

//.........这里部分代码省略.........

void determinant(ZZ& rres, const mat_ZZ& a, long deterministic)
{
    long n = a.NumRows();
    if (a.NumCols() != n)
        Error("determinant: nonsquare matrix");

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

    zz_pBak zbak;
    zbak.save();

    ZZ_pBak Zbak;
    Zbak.save();

    long instable = 1;

    long gp_cnt = 0;

    long bound = 2+DetBound(a);

    ZZ res, prod;

    clear(res);
    set(prod);


    long i;
    for (i = 0; ; i++) {
        if (NumBits(prod) > bound)
            break;

        if (!deterministic &&
                !instable && bound > 1000 && NumBits(prod) < 0.25*bound) {
            ZZ P;


            long plen = 90 + NumBits(max(bound, NumBits(res)));
            GenPrime(P, plen, 90 + 2*NumBits(gp_cnt++));

            ZZ_p::init(P);

            mat_ZZ_p A;
            conv(A, a);

            ZZ_p t;
            determinant(t, A);

            if (CRT(res, prod, rep(t), P))
                instable = 1;
            else
                break;
        }


        zz_p::FFTInit(i);
        long p = zz_p::modulus();

        mat_zz_p A;
        conv(A, a);

        zz_p t;
        determinant(t, A);

        instable = CRT(res, prod, rep(t), p);
    }

    rres = res;

    zbak.restore();
    Zbak.restore();
}

void solve1(ZZ& d_out, vec_ZZ& x_out, const mat_ZZ& A, const vec_ZZ& b)
{
   long n = A.NumRows();

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

   if (b.length() != n)
      LogicError("solve1: dimension mismatch");

   if (n == 0) {
      set(d_out);
      x_out.SetLength(0);
      return;
   }

   ZZ num_bound, den_bound;

   hadamard(num_bound, den_bound, A, b);

   if (den_bound == 0) {
      clear(d_out);
      return;
   }

   zz_pBak zbak;
   zbak.save();

   long i;
   long j;

   ZZ prod;
   prod = 1;

   mat_zz_p B;


   for (i = 0; ; i++) {
      zz_p::FFTInit(i);

      mat_zz_p AA, BB;
      zz_p dd;

      conv(AA, A);
      inv(dd, BB, AA);

      if (dd != 0) {
         transpose(B, BB);
         break;
      }

      mul(prod, prod, zz_p::modulus());
      
      if (prod > den_bound) {
         d_out = 0;
         return;
      }
   }

   long max_A_len = MaxBits(A);

   long use_double_mul1 = 0;
   long use_double_mul2 = 0;
   long double_limit = 0;

   if (max_A_len + NTL_SP_NBITS + NumBits(n) <= NTL_DOUBLE_PRECISION-1)
      use_double_mul1 = 1;

   if (!use_double_mul1 && max_A_len+NTL_SP_NBITS+2 <= NTL_DOUBLE_PRECISION-1) {
      use_double_mul2 = 1;
      double_limit = (1L << (NTL_DOUBLE_PRECISION-1-max_A_len-NTL_SP_NBITS));
   }

   long use_long_mul1 = 0;
   long use_long_mul2 = 0;
   long long_limit = 0;

   if (max_A_len + NTL_SP_NBITS + NumBits(n) <= NTL_BITS_PER_LONG-1)
      use_long_mul1 = 1;

   if (!use_long_mul1 && max_A_len+NTL_SP_NBITS+2 <= NTL_BITS_PER_LONG-1) {
      use_long_mul2 = 1;
      long_limit = (1L << (NTL_BITS_PER_LONG-1-max_A_len-NTL_SP_NBITS));
   }



   if (use_double_mul1 && use_long_mul1)
      use_long_mul1 = 0;
   else if (use_double_mul1 && use_long_mul2)
      use_long_mul2 = 0;
   else if (use_double_mul2 && use_long_mul1)
      use_double_mul2 = 0;
   else if (use_double_mul2 && use_long_mul2) {
      if (long_limit > double_limit)
         use_double_mul2 = 0;
      else
         use_long_mul2 = 0;
   }

//.........这里部分代码省略.........

void CharPoly(ZZX& gg, const mat_ZZ& a, long deterministic)
{
   long n = a.NumRows();
   if (a.NumCols() != n)
      LogicError("CharPoly: nonsquare matrix");

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


   if (n == 1) {
      ZZ t;
      SetX(gg);
      negate(t, a(1, 1));
      SetCoeff(gg, 0, t);
      return;
   }

   long bound = 2 + CharPolyBound(a);

   zz_pBak bak;
   bak.save();

   ZZ_pBak bak1;
   bak1.save();

   ZZX g;
   ZZ prod;

   clear(g);
   set(prod);

   long i;

   long instable = 1;

   long gp_cnt = 0;

   for (i = 0; ; i++) {
      if (NumBits(prod) > bound)
         break;

      if (!deterministic &&
          !instable && bound > 1000 && NumBits(prod) < 0.25*bound) {
         long plen = 90 + NumBits(max(bound, MaxBits(g)));

         ZZ P;

         GenPrime(P, plen, 90 + 2*NumBits(gp_cnt++));

         ZZ_p::init(P);
         mat_ZZ_p A;
         ZZ_pX G;
         conv(A, a);
         CharPoly(G, A);

         if (CRT(g, prod, G))
            instable = 1;
         else
            break;
      }

      zz_p::FFTInit(i);

      mat_zz_p A;
      zz_pX G;
      conv(A, a);
      CharPoly(G, A);
      instable = CRT(g, prod, G);
   }

   gg = g;

   bak.restore();
   bak1.restore();
}

static
long LLL_XD(mat_ZZ& B, mat_ZZ* U, xdouble delta, long deep, 
           LLLCheckFct check)
{
   long m = B.NumRows();
   long n = B.NumCols();

   long i, j;
   long new_m, dep, quit;
   xdouble s;
   ZZ MU;
   xdouble mu1;

   xdouble t1;
   ZZ T1;

   init_red_fudge();

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

   xdouble **B1;  // approximates B

   typedef xdouble *xdoubleptr;

   B1 = NTL_NEW_OP xdoubleptr[m+1];
   if (!B1) Error("LLL_XD: out of memory");

   for (i = 1; i <= m; i++) {
      B1[i] = NTL_NEW_OP xdouble[n+1];
      if (!B1[i]) Error("LLL_XD: out of memory");
   }

   xdouble **mu;
   mu = NTL_NEW_OP xdoubleptr[m+1];
   if (!mu) Error("LLL_XD: out of memory");

   for (i = 1; i <= m; i++) {
      mu[i] = NTL_NEW_OP xdouble[m+1];
      if (!mu[i]) Error("LLL_XD: out of memory");
   }

   xdouble *c; // squared lengths of Gramm-Schmidt basis vectors

   c = NTL_NEW_OP xdouble[m+1];
   if (!c) Error("LLL_XD: out of memory");

   xdouble *b; // squared lengths of basis vectors

   b = NTL_NEW_OP xdouble[m+1];
   if (!b) Error("LLL_XD: out of memory");



   for (i = 1; i <=m; i++)
      for (j = 1; j <= n; j++) 
         conv(B1[i][j], B(i, j));


         
   for (i = 1; i <= m; i++) {
      b[i] = InnerProduct(B1[i], B1[i], n);
   }


   new_m = ll_LLL_XD(B, U, delta, deep, check, B1, mu, b, c, m, 1, quit);
   dep = m - new_m;
   m = new_m;

   if (dep > 0) {
      // for consistency, we move all of the zero rows to the front

      for (i = 0; i < m; i++) {
         swap(B(m+dep-i), B(m-i));
         if (U) swap((*U)(m+dep-i), (*U)(m-i));
      }
   }


   // clean-up

   for (i = 1; i <= m+dep; i++) {
      delete [] B1[i];
   }

   delete [] B1;

   for (i = 1; i <= m+dep; i++) {
      delete [] mu[i];
   }

   delete [] mu;

   delete [] c;

   delete [] b;

   return m;
}

static
long BKZ_XD(mat_ZZ& BB, mat_ZZ* UU, xdouble delta, 
         long beta, long prune, LLLCheckFct check)
{
   long m = BB.NumRows();
   long n = BB.NumCols();
   long m_orig = m;
   
   long i, j;
   ZZ MU;

   xdouble t1;
   ZZ T1;
   xdouble *tp;

   init_red_fudge();

   mat_ZZ B;
   B = BB;

   B.SetDims(m+1, n);


   xdouble **B1;  // approximates B

   typedef xdouble *xdoubleptr;

   B1 = NTL_NEW_OP xdoubleptr[m+2];
   if (!B1) Error("BKZ_XD: out of memory");

   for (i = 1; i <= m+1; i++) {
      B1[i] = NTL_NEW_OP xdouble[n+1];
      if (!B1[i]) Error("BKZ_XD: out of memory");
   }

   xdouble **mu;
   mu = NTL_NEW_OP xdoubleptr[m+2];
   if (!mu) Error("BKZ_XD: out of memory");

   for (i = 1; i <= m+1; i++) {
      mu[i] = NTL_NEW_OP xdouble[m+1];
      if (!mu[i]) Error("BKZ_XD: out of memory");
   }

   xdouble *c; // squared lengths of Gramm-Schmidt basis vectors

   c = NTL_NEW_OP xdouble[m+2];
   if (!c) Error("BKZ_XD: out of memory");

   xdouble *b; // squared lengths of basis vectors

   b = NTL_NEW_OP xdouble[m+2];
   if (!b) Error("BKZ_XD: out of memory");

   xdouble cbar;

   xdouble *ctilda;
   ctilda = NTL_NEW_OP xdouble[m+2];
   if (!ctilda) Error("BKZ_XD: out of memory");

   xdouble *vvec;
   vvec = NTL_NEW_OP xdouble[m+2];
   if (!vvec) Error("BKZ_XD: out of memory");

   xdouble *yvec;
   yvec = NTL_NEW_OP xdouble[m+2];
   if (!yvec) Error("BKZ_XD: out of memory");

   xdouble *uvec;
   uvec = NTL_NEW_OP xdouble[m+2];
   if (!uvec) Error("BKZ_XD: out of memory");

   xdouble *utildavec;
   utildavec = NTL_NEW_OP xdouble[m+2];
   if (!utildavec) Error("BKZ_XD: out of memory");


   long *Deltavec;
   Deltavec = NTL_NEW_OP long[m+2];
   if (!Deltavec) Error("BKZ_XD: out of memory");

   long *deltavec;
   deltavec = NTL_NEW_OP long[m+2];
   if (!deltavec) Error("BKZ_XD: out of memory");

   mat_ZZ Ulocal;
   mat_ZZ *U;

   if (UU) {
      Ulocal.SetDims(m+1, m);
      for (i = 1; i <= m; i++)
         conv(Ulocal(i, i), 1);
      U = &Ulocal;
   }
   else
      U = 0;

   long quit;
   long new_m;
   long z, jj, kk;
//.........这里部分代码省略.........

static
long G_LLL_XD(mat_ZZ& B, mat_ZZ* U, xdouble delta, long deep, 
           LLLCheckFct check)
{
   long m = B.NumRows();
   long n = B.NumCols();

   long i, j;
   long new_m, dep, quit;
   xdouble s;
   ZZ MU;
   xdouble mu1;

   xdouble t1;
   ZZ T1;

   init_red_fudge();

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

   xdouble **B1;  // approximates B

   typedef xdouble *xdoubleptr;

   B1 = NTL_NEW_OP xdoubleptr[m+1];
   if (!B1) Error("G_LLL_XD: out of memory");

   for (i = 1; i <= m; i++) {
      B1[i] = NTL_NEW_OP xdouble[n+1];
      if (!B1[i]) Error("G_LLL_XD: out of memory");
   }

   xdouble **mu;
   mu = NTL_NEW_OP xdoubleptr[m+1];
   if (!mu) Error("G_LLL_XD: out of memory");

   for (i = 1; i <= m; i++) {
      mu[i] = NTL_NEW_OP xdouble[n+2];
      if (!mu[i]) Error("G_LLL_XD: out of memory");
   }

   xdouble **aux;
   aux = NTL_NEW_OP xdoubleptr[m+1];
   if (!aux) Error("G_LLL_XD: out of memory");

   for (i = 1; i <= m; i++) {
      aux[i] = NTL_NEW_OP xdouble[n+1];
      if (!aux[i]) Error("G_LLL_XD: out of memory");
   }

   for (i = 1; i <=m; i++)
      for (j = 1; j <= n; j++) 
         conv(B1[i][j], B(i, j));

   GivensCache_XD cache(m, n);

   new_m =
      ll_G_LLL_XD(B, U, delta, deep, check, B1, mu, aux, m, 1, quit, cache);

   dep = m - new_m;
   m = new_m;

   if (dep > 0) {
      // for consistency, we move all of the zero rows to the front

      for (i = 0; i < m; i++) {
         swap(B(m+dep-i), B(m-i));
         if (U) swap((*U)(m+dep-i), (*U)(m-i));
      }
   }


   // clean-up

   for (i = 1; i <= m+dep; i++) {
      delete [] B1[i];
   }

   delete [] B1;

   for (i = 1; i <= m+dep; i++) {
      delete [] mu[i];
   }

   delete [] mu;

   for (i = 1; i <= m+dep; i++) {
      delete [] aux[i];
   }

   delete [] aux;

   return m;
}

long LatticeSolve(vec_ZZ& x, const mat_ZZ& A, const vec_ZZ& y, long reduce)
{
   long n = A.NumRows();
   long m = A.NumCols();

   if (y.length() != m)
      Error("LatticeSolve: dimension mismatch");

   if (reduce < 0 || reduce > 2)
      Error("LatticeSolve: bad reduce parameter");

   if (IsZero(y)) {
      x.SetLength(n);
      clear(x);
      return 1;
   }

   mat_ZZ A1, U1;
   ZZ det2;
   long im_rank, ker_rank;

   A1 = A;

   im_rank = image(det2, A1, U1);
   ker_rank = n - im_rank;

   mat_ZZ A2, U2;
   long new_rank;
   long i;

   A2.SetDims(im_rank + 1, m);
   for (i = 1; i <= im_rank; i++)
      A2(i) = A1(ker_rank + i);

   A2(im_rank + 1) = y;

   new_rank = image(det2, A2, U2);

   if (new_rank != im_rank ||
      (U2(1)(im_rank+1) != 1  && U2(1)(im_rank+1) != -1))
      return 0;

   vec_ZZ x1;
   x1.SetLength(im_rank);

   for (i = 1; i <= im_rank; i++)
      x1(i) = U2(1)(i);

   if (U2(1)(im_rank+1) == 1)
      negate(x1, x1);

   vec_ZZ x2, tmp;
   x2.SetLength(n);
   clear(x2);
   tmp.SetLength(n);

   for (i = 1; i <= im_rank; i++) {
      mul(tmp, U1(ker_rank+i), x1(i));
      add(x2, x2, tmp);
   }

   if (reduce == 0) {
      x = x2;
      return 1;
   }
   else if (reduce == 1) {
      U1.SetDims(ker_rank+1, n);
      U1(ker_rank+1) = x2;
      image(det2, U1);

      x = U1(ker_rank + 1);
      return 1;
   }
   else if (reduce == 2) {
      U1.SetDims(ker_rank, n);
      LLL(det2, U1);
      U1.SetDims(ker_rank+1, n);
      U1(ker_rank+1) = x2;
      image(det2, U1);

      x = U1(ker_rank + 1);
      return 1;
   }

   return 0;
}

static
long LLL(vec_ZZ& D, mat_ZZ& B, mat_ZZ* U, long a, long b, long verbose)
{
   long m = B.NumRows();
   long n = B.NumCols();

   long force_reduce = 1;

   vec_long P;
   P.SetLength(m);

   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 ||
           SwapTest(D[P(k)], D[P(k)-1], D[P(k)-2], lam(k)(P(k)-1), a, b))) {
         force_reduce = swap(k, B, P, D, lam, U, max_k, verbose);
         k--;
      }
      else {
         force_reduce = 1;
         k++;
      }
   }

   D.SetLength(s+1);
   return s;
}

static
long LLL_XD(mat_ZZ& B, mat_ZZ* U, xdouble delta, long deep, 
           LLLCheckFct check)
{
   long m = B.NumRows();
   long n = B.NumCols();

   long i, j;
   long new_m, dep, quit;
   xdouble s;
   ZZ MU;
   xdouble mu1;

   xdouble t1;
   ZZ T1;

   init_red_fudge();

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

   Unique2DArray<xdouble> B1_store;
   B1_store.SetDimsFrom1(m+1, n+1);
   xdouble **B1 = B1_store.get();  // approximates B


   Unique2DArray<xdouble> mu_store;
   mu_store.SetDimsFrom1(m+1, m+1);
   xdouble **mu = mu_store.get();

   UniqueArray<xdouble> c_store;
   c_store.SetLength(m+1);
   xdouble *c = c_store.get(); // squared lengths of Gramm-Schmidt basis vectors

   UniqueArray<xdouble> b_store;
   b_store.SetLength(m+1);
   xdouble *b = b_store.get(); // squared lengths of basis vectors

   for (i = 1; i <=m; i++)
      for (j = 1; j <= n; j++) 
         conv(B1[i][j], B(i, j));


         
   for (i = 1; i <= m; i++) {
      b[i] = InnerProduct(B1[i], B1[i], n);
   }


   new_m = ll_LLL_XD(B, U, delta, deep, check, B1, mu, b, c, m, 1, quit);
   dep = m - new_m;
   m = new_m;

   if (dep > 0) {
      // for consistency, we move all of the zero rows to the front

      for (i = 0; i < m; i++) {
         swap(B(m+dep-i), B(m-i));
         if (U) swap((*U)(m+dep-i), (*U)(m-i));
      }
   }

   return m;
}

RR ReductionQualityChecker::computeHermiteFactor(mat_ZZ &mat) {
    RR numerator   = findShortestNormVector(mat);
    RR denominator = pow(abs(determinant(mat)), (1./mat.NumRows()));
    RR result =  numerator/denominator;
    return result;
}