C++ ZZX类代码示例

void sub(ZZX& x, const ZZ& b, const ZZX& a)
{
   long n = a.rep.length();
   if (n == 0) {
      conv(x, b);
   }
   else if (x.rep.MaxLength() == 0) {
      negate(x, a);
      add(x.rep[0], a.rep[0], b);
      x.normalize();
   }
   else {
      // ugly...b could alias a coeff of x

      ZZ *xp = x.rep.elts();
      sub(xp[0], b, a.rep[0]);
      x.rep.SetLength(n);
      xp = x.rep.elts();
      const ZZ *ap = a.rep.elts();
      long i;
      for (i = 1; i < n; i++)
         negate(xp[i], ap[i]);
      x.normalize();
   }
}

void sampleGaussian(ZZX &poly, long n, double stdev)
{
  static double const Pi=4.0*atan(1.0); // Pi=3.1415..
  static long const bignum = 0xfffffff;
  // THREADS: C++11 guarantees these are initialized only once

  if (n<=0) n=deg(poly)+1; if (n<=0) return;
  poly.SetMaxLength(n); // allocate space for degree-(n-1) polynomial
  for (long i=0; i<n; i++) SetCoeff(poly, i, ZZ::zero());

  // Uses the Box-Muller method to get two Normal(0,stdev^2) variables
  for (long i=0; i<n; i+=2) {
    double r1 = (1+RandomBnd(bignum))/((double)bignum+1);
    double r2 = (1+RandomBnd(bignum))/((double)bignum+1);
    double theta=2*Pi*r1;
    double rr= sqrt(-2.0*log(r2))*stdev;

    assert(rr < 8*stdev); // sanity-check, no more than 8 standard deviations

    // Generate two Gaussians RV's, rounded to integers
    long x = (long) floor(rr*cos(theta) +0.5);
    SetCoeff(poly, i, x);
    if (i+1 < n) {
      x = (long) floor(rr*sin(theta) +0.5);
      SetCoeff(poly, i+1, x);
    }
  }
  poly.normalize(); // need to call this after we work on the coeffs
}

ZZX RandPoly(long n,const ZZ& p)
{ 
  ZZX F; F.SetMaxLength(n);
  ZZ p2;  p2=p>>1;
  for (long i=0; i<n; i++)
    { SetCoeff(F,i,RandomBnd(p)-p2); }
  return F;
}

void sampleSmall(ZZX &poly, long n)
{
  if (n<=0) n=deg(poly)+1; if (n<=0) return;
  poly.SetMaxLength(n); // allocate space for degree-(n-1) polynomial

  for (long i=0; i<n; i++) {    // Chosse coefficients, one by one
    long u = lrand48();
    if (u&1) {                 // with prob. 1/2 choose between -1 and +1
      u = (u & 2) -1;
      SetCoeff(poly, i, u);
    }
    else SetCoeff(poly, i, 0); // with ptob. 1/2 set to 0
  }
  poly.normalize(); // need to call this after we work on the coeffs
}

void SingleCRT::toPoly(ZZX& poly, const IndexSet& s) const
{
  IndexSet s1 = map.getIndexSet() & s;

  if (card(s1) == 0) {
    clear(poly);
    return;
  }

  ZZ p = to_ZZ(context.ithPrime(s1.first()));  // the first modulus

  poly = map[s1.first()];  // Get poly modulo the first prime

  vec_ZZ& vp = poly.rep;

  // ensure that coeficient vector is of size phi(m) with entries in [-p/2,p/2]
  long phim = context.zMstar.phiM();
  long vpLength = vp.length();
  if (vpLength<phim) { // just in case of leading zeros in poly
    vp.SetLength(phim);
    for (long j=vpLength; j<phim; j++) vp[j]=0;
  }
  ZZ p_over_2 = p/2;
  for (long j=0; j<phim; j++) if (vp[j] > p_over_2) vp[j] -= p;

  // do incremental integer CRT for other levels  
  for (long i = s1.next(s1.first()); i <= s1.last(); i = s1.next(i)) {
    long q = context.ithPrime(i);       // the next modulus

    // CRT the coefficient vectors of poly and current
    intVecCRT(vp, p, map[i].rep, q);    // defined in the module NumbTh
    p *= q;     // update the modulus
  }
  poly.normalize(); // need to call this after we work on the coeffs
}

void PlainMul(ZZX& x, const ZZX& a, const ZZX& b)
{
   if (&a == &b) {
      PlainSqr(x, a);
      return;
   }

   long da = deg(a);
   long db = deg(b);

   if (da < 0 || db < 0) {
      clear(x);
      return;
   }

   long d = da+db;



   const ZZ *ap, *bp;
   ZZ *xp;

   ZZX la, lb;

   if (&x == &a) {
      la = a;
      ap = la.rep.elts();
   }
   else
      ap = a.rep.elts();

   if (&x == &b) {
      lb = b;
      bp = lb.rep.elts();
   }
   else
      bp = b.rep.elts();

   x.rep.SetLength(d+1);

   xp = x.rep.elts();

   long i, j, jmin, jmax;
   ZZ t, accum;

   for (i = 0; i <= d; i++) {
      jmin = max(0, i-db);
      jmax = min(da, i);
      clear(accum);
      for (j = jmin; j <= jmax; j++) {
	 mul(t, ap[j], bp[i-j]);
	 add(accum, accum, t);
      }
      xp[i] = accum;
   }
   x.normalize();
}

void KarSqr(ZZX& c, const ZZX& a)
{
   if (IsZero(a)) {
      clear(c);
      return;
   }

   vec_ZZ mem;

   const ZZ *ap;
   ZZ *cp;

   long sa = a.rep.length();

   if (&a == &c) {
      mem = a.rep;
      ap = mem.elts();
   }
   else
      ap = a.rep.elts();

   c.rep.SetLength(sa+sa-1);
   cp = c.rep.elts();

   long maxa, xover;

   maxa = MaxBits(a);

   xover = 2;

   if (sa < xover)
      PlainSqr(cp, ap, sa);
   else {
      /* karatsuba */

      long n, hn, sp, depth;

      n = sa;
      sp = 0;
      depth = 0;
      do {
         hn = (n+1) >> 1;
         sp += hn+hn+hn - 1;
         n = hn;
         depth++;
      } while (n >= xover);

      ZZVec stk;
      stk.SetSize(sp,
         ((2*maxa + NumBits(sa) + 2*depth + 10)
          + NTL_ZZ_NBITS-1)/NTL_ZZ_NBITS);

      KarSqr(cp, ap, sa, stk.elts());
   }

   c.normalize();
}

void add(ZZX& x, const ZZX& a, long b)
{
   if (a.rep.length() == 0) {
      conv(x, b);
   }
   else {
      if (&x != &a) x = a;
      add(x.rep[0], x.rep[0], b);
      x.normalize();
   }
}

void PlainSqr(ZZX& x, const ZZX& a)
{
   long da = deg(a);

   if (da < 0) {
      clear(x);
      return;
   }

   long d = 2*da;

   const ZZ *ap;
   ZZ *xp;

   ZZX la;

   if (&x == &a) {
      la = a;
      ap = la.rep.elts();
   }
   else
      ap = a.rep.elts();


   x.rep.SetLength(d+1);

   xp = x.rep.elts();

   long i, j, jmin, jmax;
   long m, m2;
   ZZ t, accum;

   for (i = 0; i <= d; i++) {
      jmin = max(0, i-da);
      jmax = min(da, i);
      m = jmax - jmin + 1;
      m2 = m >> 1;
      jmax = jmin + m2 - 1;
      clear(accum);
      for (j = jmin; j <= jmax; j++) {
	 mul(t, ap[j], ap[i-j]);
	 add(accum, accum, t);
      }
      add(accum, accum, accum);
      if (m & 1) {
	 sqr(t, ap[jmax + 1]);
	 add(accum, accum, t);
      }

      xp[i] = accum;
   }

   x.normalize();
}

void sampleHWt(ZZX &poly, long Hwt, long n)
{
  if (n<=0) n=deg(poly)+1; if (n<=0) return;
  clear(poly);          // initialize to zero
  poly.SetMaxLength(n); // allocate space for degree-(n-1) polynomial

  long b,u,i=0;
  if (Hwt>n) Hwt=n;
  while (i<Hwt) {  // continue until exactly Hwt nonzero coefficients
    u=lrand48()%n; // The next coefficient to choose
    if (IsZero(coeff(poly,u))) { // if we didn't choose it already
      b = lrand48()&2; // b random in {0,2}
      b--;             //   random in {-1,1}
      SetCoeff(poly,u,b);

      i++; // count another nonzero coefficient
    }
  }
  poly.normalize(); // need to call this after we work on the coeffs
}

void sampleUniform(ZZX& poly, const ZZ& B, long n)
{
  if (n<=0) n=deg(poly)+1; if (n<=0) return;
  if (B <= 0) {
    clear(poly);
    return;
  }

  poly.SetMaxLength(n); // allocate space for degree-(n-1) polynomial

  ZZ UB, tmp;

  UB =  2*B + 1;
  for (long i = 0; i < n; i++) {
    RandomBnd(tmp, UB);
    tmp -= B; 
    poly.rep[i] = tmp;
  }

  poly.normalize();
}

// multiply the polynomial f by the integer a modulo q
void MulMod(ZZX& out, const ZZX& f, long a, long q, bool abs/*default=true*/)
{
  // ensure that out has the same degree as f
  out.SetMaxLength(deg(f)+1);               // allocate space if needed
  if (deg(out)>deg(f)) trunc(out,out,deg(f)+1); // remove high degrees

  for (int i=0; i<=deg(f); i++) { 
    int c = rem(coeff(f,i), q);
    c = MulMod(c, a, q); // returns c \in [0,q-1]
    if (!abs && c >= q/2)
      c -= q;
    SetCoeff(out,i,c);
  }
}

void SetCoeff(ZZX& x, long i)
{
   long j, m;

   if (i < 0)
      Error("coefficient index out of range");

   if (NTL_OVERFLOW(i, 1, 0))
      Error("overflow in SetCoeff");

   m = deg(x);

   if (i > m) {
      x.rep.SetLength(i+1);
      for (j = m+1; j < i; j++)
         clear(x.rep[j]);
   }
   set(x.rep[i]);
   x.normalize();
}

void PolyRed(ZZX& out, const ZZX& in, long q, bool abs)
{
  // ensure that out has the same degree as in
  out.SetMaxLength(deg(in)+1);               // allocate space if needed
  if (deg(out)>deg(in)) trunc(out,out,deg(in)+1); // remove high degrees

  long q2; q2=q>>1;
  for (long i=0; i<=deg(in); i++)
    { long c=coeff(in,i)%q;
      if (abs)
        { if (c<0) { c=c+q; } }
      else if (q==2)
        { if (coeff(in,i)<0) { c=-c; } }
      else
        { if (c>=q2)  { c=c-q; }
          else if (c<-q2) { c=c+q; }
	}
      SetCoeff(out,i,c);
    }
}

//FIXME: both the reduction from powerful to the individual primes and
//  the CRT back to poly can be made more efficient
void PowerfulDCRT::powerfulToZZX(ZZX& poly, const Vec<ZZ>& powerful,
				 IndexSet set) const
{
  zz_pBak bak; bak.save(); // backup NTL's current modulus

  if (empty(set)) set = IndexSet(0, pConvVec.length()-1);

  clear(poly);
  //  poly.SetLength(powerful.length());
  ZZ product = conv<ZZ>(1L);
  for (long i = set.first(); i <= set.last(); i = set.next(i)) {
    pConvVec[i].restoreModulus();
    //    long newPrime = zz_p::modulus();

    HyperCube<zz_p> oneRowPwrfl(indexes.shortSig);
    conv(oneRowPwrfl.getData(), powerful); // reduce and convert to Vec<zz_p>

    zz_pX oneRowPoly;
    pConvVec[i].powerfulToPoly(oneRowPoly, oneRowPwrfl);
    CRT(poly, product, oneRowPoly);                   // NTL :-)
  }
  poly.normalize();
}