C++ vec_ZZ::SetLength方法代码示例

void negate(vec_ZZ& x, const vec_ZZ& a)
{
   long n = a.length();
   x.SetLength(n);
   long i;
   for (i = 0; i < n; i++)
      negate(x[i], a[i]);
}

void conv(vec_ZZ& x, const vec_zz_p& a)
{
   long n = a.length();
   x.SetLength(n);
   long i;
   for (i = 0; i < n; i++)
      x[i] = rep(a[i]);
}

void mul(vec_ZZ& x, const vec_ZZ& a, long b)
{
   long n = a.length();
   x.SetLength(n);
   long i;
   for (i = 0; i < n; i++)
      mul(x[i], a[i], b);
}

void sub(vec_ZZ& x, const vec_ZZ& a, const vec_ZZ& b)
{
   long n = a.length();
   if (b.length() != n) LogicError("vector sub: dimension mismatch");
   x.SetLength(n);
   long i;
   for (i = 0; i < n; i++)
      sub(x[i], a[i], b[i]);
}

static void MulSubDiv(vec_ZZ& c, const vec_ZZ& c1, const vec_ZZ& c2,
                      const ZZ& x, const ZZ& y, const ZZ& z)

// c = (x*c1 + y*c2)/z

{
   long n = c1.length();
   if (c2.length() != n) Error("MulSubDiv: length mismatch");
   c.SetLength(n);

   long i;
   for (i = 1; i <= n; i++)
      MulSubDiv(c(i), c1(i), c2(i), x, y, z);
}

void VectorCopy(vec_ZZ& x, const vec_ZZ& a, long n)
{
   if (n < 0) LogicError("VectorCopy: negative length");
   if (NTL_OVERFLOW(n, 1, 0)) ResourceError("overflow in VectorCopy");

   long m = min(n, a.length());

   x.SetLength(n);

   long i;

   for (i = 0; i < m; i++)
      x[i] = a[i];

   for (i = m; i < n; i++)
      clear(x[i]);
}

//Ideal lattice challenge problem generator 
//modified from generate_ideal.cpp in Ideal lattice Challenge web page
void gen_idealsvpchallenge(mat_ZZ& B,int index,ZZ seed,vec_ZZ& phivec) {

    ZZX phi=find_cyclotomic(index);  
    long n=deg(phi);
    ZZ det=find_determinant(index,10*n,seed);
    ZZ alpha=find_unity_root(index,det,phi);
    B.SetDims(n,n); 
    clear(B);
    B(1,1) = det;
    for (long i=2; i<=n; i++)
      {
	B(i,1)=det-PowerMod(alpha,i-1,det);
	B(i,i)=1;
      }
    
    phivec.SetLength(n+1);
    for (int i=0;i<=n;i++) phivec[i] = coeff(phi,i);
}

static
void MixedMul(vec_ZZ& x, const vec_zz_p& a, const mat_ZZ& B)
{
    long n = B.NumRows();
    long l = B.NumCols();

    if (n != a.length())
        Error("matrix mul: dimension mismatch");

    x.SetLength(l);

    long i, k;
    ZZ acc, tmp;

    for (i = 1; i <= l; i++) {
        clear(acc);
        for (k = 1; k <= n; k++) {
            mul(tmp, B(k, i), rep(a(k)));
            add(acc, acc, tmp);
        }
        x(i) = acc;
    }
}

static
void mul_aux(vec_ZZ& x, const mat_ZZ& A, const vec_ZZ& b)
{
    long n = A.NumRows();
    long l = A.NumCols();

    if (l != b.length())
        Error("matrix mul: dimension mismatch");

    x.SetLength(n);

    long i, k;
    ZZ acc, tmp;

    for (i = 1; i <= n; i++) {
        clear(acc);
        for (k = 1; k <= l; k++) {
            mul(tmp, A(i,k), b(k));
            add(acc, acc, tmp);
        }
        x(i) = acc;
    }
}

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

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

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

//.........这里部分代码省略.........
    std::cout<<"m="<<m<<std::endl;
    std::cout<<"q="<<q<<std::endl;
    // step 6
    EC_pCurve curve;
    ZZ_p A,B;
    find_curve(A,B,to_long(D),N,HCP);
    SetCoeff(curve,1,A);
    SetCoeff(curve,0,B);
    std::cout<<"curve="<<curve<<std::endl;
    // step 7
    ZZ_p g;
    do {
      random(g);
      if (Jacobi(rep(g),N)!=-1)
	continue;
      if ((N%3)==1) {
	// always happens in the case D==-3
	ZZ_p t;
	power(t,g,(N-1)/3);
	if (IsOne(t))
	  continue;
	//std::cout<<"t="<<t<<std::endl;
      }
      break;
    } while (true);
    //std::cout<<"g="<<g<<std::endl;
    long k=0;
    EC_p::init(curve);
    EC_p P,P1,P2;
    do {
      // step 8
      random(P);
      if (!IsValid(P)) 
	return false;
      // step 9
      //std::cout<<"P="<<P<<std::endl;
      mul(P2,P,m/q);
      //std::cout<<"P2="<<P2<<std::endl;
      // NOTE: if IsZero(P2) just choose a new point, not a new curve!!!
      if (!IsZero(P2)) {
	mul(P1,P2,q);
	//std::cout<<"P1="<<P1<<std::endl;
	if (IsZero(P1)) 
	  break;
      }
      // step 10
      ++k;
      if (D==-3) {
	if (k>=6) {
	  std::cout<<"P="<<P<<std::endl;
	  std::cout<<"P1="<<P1<<std::endl;
	  std::cout<<"P2="<<P2<<std::endl;
	  return false;
	}
	SetCoeff(curve,0,coeff(curve,0)*g);
      }
      else if (D==-4) {
	if (k>=4) {
	  std::cout<<"P="<<P<<std::endl;
	  std::cout<<"P1="<<P1<<std::endl;
	  std::cout<<"P2="<<P2<<std::endl;
	  return false;
	}
	SetCoeff(curve,1,coeff(curve,1)*g);
      }
      else {
	if (k>=2) {
	  std::cout<<"P="<<P<<std::endl;
	  std::cout<<"P1="<<P1<<std::endl;
	  std::cout<<"P2="<<P2<<std::endl;
	  return false;
	}
	SetCoeff(curve,1,coeff(curve,1)*sqr(g));
	SetCoeff(curve,0,coeff(curve,0)*sqr(g)*g);
      }
      EC_p::init(curve);
    } while (true);

    // check that P2 has an affine representation
    ZZ G;
    GCD(G,rep(P2.Z),N);
    if (!IsOne(G))
      return false;

    // step 13
    ZZ_p Px,Py;
    P.affine(Px,Py);
    cert.SetLength(7);
    cert[0] = N;
    cert[1] = q;
    cert[2] = m/q;
    cert[3] = coeff(curve,1)==-1 ? to_ZZ(-1) : rep(coeff(curve,1));
    cert[4] = coeff(curve,0)==-1 ? to_ZZ(-1) : rep(coeff(curve,0));
    cert[5] = rep(Px);
    cert[6] = rep(Py);
    return true;
  }
  Error("ProvePrime: ran out of discriminants");
  return false;
}