C++ SetCoeff函数代码示例

// choose random curve
// if IsZero(Z), then X is non-trivial factor of ZZ_p::modulus()
void ECM_random_curve(EC_pCurve& curve, ZZ_p& X, ZZ_p& Z) {
  ZZ sigma;
  RandomBnd(sigma,ZZ_p::modulus()-6);
  sigma+=6;
  ZZ_p u,v;
  u = sqr(to_ZZ_p(sigma))-5;
  v = 4*to_ZZ_p(sigma);
  ZZ_p C,Cd;
  C = (v-u)*sqr(v-u)*(3*u+v);
  Cd = 4*u*sqr(u)*v;
  // make sure Cd is invertible
  ZZ Cinv;
  if (InvModStatus(Cinv,rep(Cd),ZZ_p::modulus())!=0) {
    conv(X,Cinv);
    clear(Z);
    return;
  }
  C*=to_ZZ_p(Cinv);
  C-=2;

  // random curve
  ZZ_pX f;
  SetCoeff(f,3);
  SetCoeff(f,2,C);
  SetCoeff(f,1);
  conv(curve,f);
  curve.SetRepresentation(curve.MONTGOMERY);

  // initial point
  mul(X,u,sqr(u));
  mul(Z,v,sqr(v));
}

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
}

divisor& divisor::random(divdeg_t dgr){ // Underlying curve is not touched
  assert(s_hcurve.is_valid_curve());

  // A random valid divisor is generated by the following algorithm:
  // generate a degree 1 divisor [x - a1, b1] by choosing a1 by random
  // then trying to solve quadratic equation
  // x^2 + h(a1)*x - f(a1) for b1.
  // Note that finding a root of an equation by calling routine
  // FindRoot(root, poly) may go into an infinite loop if poly does
  // not split completely.  We avoid this by calling irreducibility 
  // test routine DetIrredTest(poly).  After a degree 1 divisor is
  // found, this divisor is doubled by calling add_cantor() to return
  // a degree 2 polynomial.

  field_t a1, b1, f_of_a1, h_of_a1;

  poly_t poly;  // polynomial x^2 + h(a1)*x - f(a1)

  SetCoeff(poly, 2); // set degree 2 leading term to 1

  do {
    do{
      NTL_NNS random(a1);

      eval(f_of_a1, s_hcurve.get_f(), a1);

      eval(h_of_a1, s_hcurve.get_h(), a1);

      SetCoeff(poly, 1, h_of_a1);
      SetCoeff(poly, 0, - f_of_a1);

    } while ( DetIrredTest(poly) );

  FindRoot(b1, poly);

    // Set upoly = x - a1
    SetX(upoly);
    SetCoeff(upoly, 0, -a1);

    // Set vpoly = b1
    vpoly = b1;

    update();
  } while (*this == -*this); // Avoid getting unit after doubling

  // return a degree one divisor if dgr = 1
  if (dgr == DEGREE_1)
    return *this;

  // Double the degree 1 divisor to get a degree 2 divisor, otherwise
  add_cantor(*this, *this, *this);

  if (is_valid_divisor())
    return *this;

  cerr << "Random divisor failed to generate" << endl;
  abort();
}

void sample_NTL_poly_div1(unsigned long length, unsigned long bits,
                          void* arg, unsigned long count)
{

       
    ZZX poly1;
    ZZX poly2;
    ZZX poly3;
    ZZ a;
    poly1.SetMaxLength(length);
    poly2.SetMaxLength(length);
    poly3.SetMaxLength(2*length-1);


   
   
   unsigned long r_count;    // how often to generate new random data
   
   if (count >= 10000) r_count = 100;
   else if (count >= 100) r_count = 10;
   else if (count >= 20) r_count = 4;
   else if (count >= 8) r_count = 2;
   else r_count = 1;
   
   unsigned long i;
   for (i = 0; i < count; i++)
   {
      if (i%r_count == 0)
      {
	    do
        {
           unsigned long j;
           for (j = 0; j < length; j++)
		   {
		      RandomBits(a,bits);
		      SetCoeff(poly1,j,a);
		   }
        } while (IsZero(poly1));
        unsigned long j;
        for (j = 0; j < length; j++)
		{
           RandomBits(a,bits);
		   SetCoeff(poly2,j,a);
		}
      }
      
      mul(poly3, poly1, poly2);
      prof_start();
      unsigned long count2;
      for (count2 = 0; count2 < r_count; count2++)
      {
         divide(poly2, poly3, poly1);
      }
      prof_stop();
      
      i += (r_count-1);
   }  
}

NTL_CLIENT

int main()
{
   zz_p::init(17);

   zz_pX P;
   BuildIrred(P, 10);

   zz_pE::init(P);

   zz_pEX f, g, h;

   random(f, 20);
   SetCoeff(f, 20);

   random(h, 20);

   g = MinPolyMod(h, f);

   if (deg(g) < 0) Error("bad zz_pEXTest (1)");
   if (CompMod(g, h, f) != 0)
      Error("bad zz_pEXTest (2)");


   
   vec_pair_zz_pEX_long v;

   long i;
   for (i = 0; i < 5; i++) {
      long n = RandomBnd(20)+1;
      cerr << n << " ";

      random(f, n);
      SetCoeff(f, n);

      v = CanZass(f);

      g = mul(v);
      if (f != g) cerr << "oops1\n";

      long i;
      for (i = 0; i < v.length(); i++)
         if (!DetIrredTest(v[i].a))
            Error("bad zz_pEXTest (3)");


   }

   cerr << "\n";

   cerr << "zz_pEXTest OK\n";
}

NTL_CLIENT


/*------------------------------------------------------------*/
/* if opt = 1, runs a check                                   */
/* else, runs timings                                         */
/*------------------------------------------------------------*/
void check(int opt){

  CTFT_multipliers mult = CTFT_multipliers(0);

  for (long i = 1; i < 770; i+=1){

    long j = 0;
    while ((1 << j) < (2*i+1)){
      CTFT_init_multipliers(mult, j);
      j++;
    }

    zz_pX A, B1, B2, C;
    for (long j = 0; j < i; j++){
      SetCoeff(A, j, random_zz_p());
      SetCoeff(B1, j, random_zz_p());
      SetCoeff(B2, j, random_zz_p());
    }
    for (long j = 0; j < 2*i-1; j++)
      SetCoeff(C, j, random_zz_p());

    if (opt == 1){
      CTFT_middle_product(B1, A, C, i, mult);
      
      B2 = A*C;
      
      for (long j = 0; j < i; j++)
	assert (coeff(B1, j) == coeff(B2, j+i-1));

      cout << i << endl;
		
    }
    else{
      cout << i;

      double u = GetTime();
      for (long j = 0; j < 10000; j++)
     	CTFT_middle_product(B1, A, C, i, mult);
      u = GetTime()-u;
      cout << " " << u;

      cout << endl;
    }
  }
}

void sample_NTL_poly_div2(unsigned long length, unsigned long bits,
                          void* arg, unsigned long count)
{

       
    ZZX poly1;
    ZZX poly2;
    ZZX poly3;
    ZZ a;
    poly1.SetMaxLength(length);
    poly2.SetMaxLength(length);
    poly3.SetMaxLength(2*length-1);


   
   
   unsigned long r_count;    // how often to generate new random data
   
   if (count >= 1000) r_count = 100;
   else if (count >= 100) r_count = 10;
   else if (count >= 20) r_count = 5;
   else if (count >= 8) r_count = 2;
   else r_count = 1;
   
   unsigned long i;
   for (i = 0; i < count; i++)
   {
      if (i%r_count == 0)
      {
	    unsigned long j;
	    for (j = 0; j<length-1; j++)
		{
		RandomBits(a,bits);
		SetCoeff(poly1,j,a);
		}
		SetCoeff(poly1,length-1,1);
	    unsigned long j;
	    for (j = 0; j<2*length-1; j++)
		{
		RandomBits(a,bits);
		SetCoeff(poly3,j,a);
		}
      }
       prof_start();
       div(poly2, poly3, poly1);
       prof_stop();
   }
   
   
}

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 EvaluationHashFunction::generateIrredPolynomial(void)
{
	// An irreducible polynomial to be used as the modulus
	GF2X p; 

	//generate the irreducible polynomial x^64+x^4+x^3+x+1 to work with 
	SetCoeff(p, 64);
	SetCoeff(p, 4);
	SetCoeff(p, 3);
	SetCoeff(p, 1);
	SetCoeff(p, 0);

	//init the field with the newly generated polynomial.
	GF2E::init(p);
}

void SetCoeff(GF2X& x, long i, long val)
{
   if (i < 0) {
      LogicError("SetCoeff: negative index");
      return; // NOTE: helps the compiler optimize
   }

   val = val & 1;

   if (val) {
      SetCoeff(x, i);
      return;
   }

   // we want to clear position i

   long n;

   n = x.xrep.length();
   long wi = i/NTL_BITS_PER_LONG;

   if (wi >= n) 
      return;

   long bi = i - wi*NTL_BITS_PER_LONG;

   x.xrep[wi] &= ~(1UL << bi);
   if (wi == n-1 && !x.xrep[wi]) x.normalize();
}

void SetCoeff(GF2X& x, long i, long val)
{
   if (i < 0) {
      Error("SetCoeff: negative index");
      return;
   }

   val = val & 1;

   if (val) {
      SetCoeff(x, i);
      return;
   }

   // we want to clear position i

   long n;

   n = x.xrep.length();
   long wi = i/newNTL_BITS_PER_LONG;

   if (wi >= n) 
      return;

   long bi = i - wi*newNTL_BITS_PER_LONG;

   x.xrep[wi] &= ~(1UL << bi);
   if (wi == n-1) x.normalize();
}

void FindRoot(GF2E& root, const GF2EX& ff)
// finds a root of ff.
// assumes that ff is monic and splits into distinct linear factors

{
   GF2EXModulus F;
   GF2EX h, h1, f;
   GF2E r;

   f = ff;
   
   if (!IsOne(LeadCoeff(f)))
      Error("FindRoot: bad args");

   if (deg(f) == 0)
      Error("FindRoot: bad args");


   while (deg(f) > 1) {
      build(F, f);
      random(r);
      clear(h);
      SetCoeff(h, 1, r);
      TraceMap(h, h, F);
      GCD(h, h, f);
      if (deg(h) > 0 && deg(h) < deg(f)) {
         if (deg(h) > deg(f)/2)
            div(f, f, h);
         else
            f = h;
      }
   }
 
   root = ConstTerm(f);
}

void *PolynomialLeader::receivePolynomials(void *ctx_tmp) {

    polynomial_rcv_ctx* ctx = (polynomial_rcv_ctx*) ctx_tmp;

    GF2E::init(*(ctx->irreduciblePolynomial));

    uint32_t coefficientSize = ctx->maskbytelen;

    uint8_t *masks = (uint8_t* )calloc(ctx->maskbytelen*ctx->numOfHashFunction*ctx->numPolynomPerHashFunc*ctx->polynomSize,sizeof(uint8_t));

    ctx->sock->Receive(masks, ctx->maskbytelen*ctx->numOfHashFunction*ctx->numPolynomPerHashFunc*ctx->polynomSize);

    for (uint32_t i = 0; i < ctx->numOfHashFunction; i++) {
        PRINT() << "receive polynomial " << i << std::endl;
        for (uint32_t k =0; k < ctx->numPolynomPerHashFunc; k++) {
            ctx->polynoms[i*ctx->numPolynomPerHashFunc+k] = new NTL::GF2EX();
            for (uint32_t j=0; j < ctx->polynomSize; j++) {

                NTL::GF2E coeffElement = PolynomialUtils::convertBytesToGF2E(
                        masks+((i*ctx->numPolynomPerHashFunc+k)*ctx->polynomSize+j)*coefficientSize, coefficientSize);
                SetCoeff(*ctx->polynoms[i*ctx->numPolynomPerHashFunc+k], j, coeffElement);
            }
        }
    }

    free(masks);
}

NTL_CLIENT

// Compute a degree-p polynomial poly(x) s.t. for any t<e and integr z of the
// form z = z0 + p^t*z1 (with 0<=z0<p), we have poly(z) = z0 (mod p^{t+1}).
//
// We get poly(x) by interpolating a degree-(p-1) polynomial poly'(x)
// s.t. poly'(z0)=z0 - z0^p (mod p^e) for all 0<=z0<p, and then setting
// poly(x) = x^p + poly'(x).
static void buildDigitPolynomial(ZZX& result, long p, long e)
{
  if (p<2 || e<=1) return; // nothing to do
  FHE_TIMER_START;
  long p2e = power_long(p,e); // the integer p^e

  // Compute x - x^p (mod p^e), for x=0,1,...,p-1
  vec_long x(INIT_SIZE, p);
  vec_long y(INIT_SIZE, p);
  long bottom = -(p/2);
  for (long j=0; j<p; j++) {
    long z = bottom+j;
    x[j] = z;
    y[j] = z-PowerMod((z < 0 ? z + p2e : z), p, p2e);  // x - x^p (mod p^e)

    while (y[j] > p2e/2)         y[j] -= p2e;
    while (y[j] < -(p2e/2))      y[j] += p2e;
  }
  interpolateMod(result, x, y, p, e);
  assert(deg(result)<p); // interpolating p points, should get deg<=p-1
  SetCoeff(result, p);   // return result = x^p + poly'(x)
  //  cerr << "# digitExt mod "<<p<<"^"<<e<<"="<<result<<endl;
  FHE_TIMER_STOP;
}

zz_pEExtraInfoT::zz_pEExtraInfoT(int precompute_inverses,
    int precompute_square_roots, int precompute_legendre_char,
    int precompute_pth_frobenius_map)
{
    int p = zz_p::modulus();
    q = to_long(zz_pE::cardinality());

    ref_count = 1;

    inv_table = precompute_inverses ? new invTable(q) : NULL;
    root_table = precompute_square_roots ? new rootTable(q) : NULL;
    legendre_table = precompute_legendre_char ? new legendreChar(q) : NULL;
    frob_map = precompute_pth_frobenius_map ? new frobeniusMap(q) : NULL;

    // precompute a non-square in F_q
    zz_pE   x, y;
    do {
	x = random_zz_pE();
    } while(x == 0 || legendre_char(x) == 1);
    non_square = x;

    // precompute image of basis 1,x^1,...,x^{d-1} under Frobenius
    frob_of_basis = new zz_pE[zz_pE::degree()];
    x = 0;
    SetCoeff(x.LoopHole(), 1, 1);
    frob_of_basis[0] = 1;
    for (int i = 1; i < zz_pE::degree(); i++) {
        power(y, x, i);
        power(frob_of_basis[i], y, p);
    }
}