zz_px1.cpp

数值算法库for Windows

void ProbMinPolyMod(ZZ_pX& h, const ZZ_pX& g, const ZZ_pXModulus& F, long m) 
{
   long n = F.n;
   if (m < 1 || m > n) Error("ProbMinPoly: bad args");
   
   long i;
   vec_ZZ_p R(INIT_SIZE, n);

   for (i = 0; i < n; i++) random(R[i]);
   DoMinPolyMod(h, g, F, m, R);
}

void MinPolyMod(ZZ_pX& hh, const ZZ_pX& g, const ZZ_pXModulus& F, long m) 
{
   ZZ_pX h, h1;
   long n = F.n;
   if (m < 1 || m > n) Error("MinPoly: bad args");

   /* probabilistically compute min-poly */

   ProbMinPolyMod(h, g, F, m);
   if (deg(h) == m) { hh = h; return; }
   CompMod(h1, h, g, F);
   if (IsZero(h1)) { hh = h; return; }

   /* not completely successful...must iterate */

   long i;

   ZZ_pX h2, h3;
   ZZ_pXMultiplier H1;
   vec_ZZ_p R(INIT_SIZE, n);

   for (;;) {
      R.SetLength(n);
      for (i = 0; i < n; i++) random(R[i]);
      build(H1, h1, F);
      UpdateMap(R, R, H1, F);
      DoMinPolyMod(h2, g, F, m-deg(h), R);

      mul(h, h, h2);
      if (deg(h) == m) { hh = h; return; }
      CompMod(h3, h2, g, F);
      MulMod(h1, h3, H1, F);
      if (IsZero(h1)) { hh = h; return; }
   }
}



void IrredPolyMod(ZZ_pX& h, const ZZ_pX& g, const ZZ_pXModulus& F, long m) 
{
   vec_ZZ_p R(INIT_SIZE, 1);
   if (m < 1 || m > F.n) Error("IrredPoly: bad args");

   set(R[0]);
   DoMinPolyMod(h, g, F, m, R);
}



void diff(ZZ_pX& x, const ZZ_pX& a)
{
   long n = deg(a);
   long i;

   if (n <= 0) {
      clear(x);
      return;
   }

   if (&x != &a)
      x.rep.SetLength(n);

   for (i = 0; i <= n-1; i++) {
      mul(x.rep[i], a.rep[i+1], i+1);
   }

   if (&x == &a)
      x.rep.SetLength(n);

   x.normalize();
}


void MakeMonic(ZZ_pX& x)
{
   if (IsZero(x))
      return;

   if (IsOne(LeadCoeff(x)))
      return;

   ZZ_p t;

   inv(t, LeadCoeff(x));
   mul(x, x, t);
}


      
void PlainMulTrunc(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b, long n)
{
   ZZ_pX y;
   mul(y, a, b);
   trunc(x, y, n);
}


void FFTMulTrunc(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b, long n)
{
   if (IsZero(a) || IsZero(b)) {
      clear(x);
      return;
   }

   long d = deg(a) + deg(b);
   if (n > d + 1)
      n = d + 1;

   long k = NextPowerOfTwo(d + 1);
   FFTRep R1(INIT_SIZE, k), R2(INIT_SIZE, k);

   ToFFTRep(R1, a, k);
   ToFFTRep(R2, b, k);
   mul(R1, R1, R2);
   FromFFTRep(x, R1, 0, n-1);
}

void MulTrunc(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b, long n)
{
   if (n < 0) Error("MulTrunc: bad args");

   if (deg(a) <= NTL_ZZ_pX_FFT_CROSSOVER || deg(b) <= NTL_ZZ_pX_FFT_CROSSOVER)
      PlainMulTrunc(x, a, b, n);
   else
      FFTMulTrunc(x, a, b, n);
}

      
void PlainSqrTrunc(ZZ_pX& x, const ZZ_pX& a, long n)
{
   ZZ_pX y;
   sqr(y, a);
   trunc(x, y, n);
}


void FFTSqrTrunc(ZZ_pX& x, const ZZ_pX& a, long n)
{
   if (IsZero(a)) {
      clear(x);
      return;
   }

   long d = 2*deg(a);
   if (n > d + 1)
      n = d + 1;

   long k = NextPowerOfTwo(d + 1);
   FFTRep R1(INIT_SIZE, k);

   ToFFTRep(R1, a, k);
   mul(R1, R1, R1);
   FromFFTRep(x, R1, 0, n-1);
}

void SqrTrunc(ZZ_pX& x, const ZZ_pX& a, long n)
{
   if (n < 0) Error("SqrTrunc: bad args");

   if (deg(a) <= NTL_ZZ_pX_FFT_CROSSOVER)
      PlainSqrTrunc(x, a, n);
   else
      FFTSqrTrunc(x, a, n);
}


void FastTraceVec(vec_ZZ_p& S, const ZZ_pX& f)
{
   long n = deg(f);

   if (n <= 0) 
      Error("FastTraceVec: bad args");

   if (n == 0) {
      S.SetLength(0);
      return;
   }

   if (n == 1) {
      S.SetLength(1);
      set(S[0]);
      return;
   }
   
   long i;
   ZZ_pX f1;

   f1.rep.SetLength(n-1);
   for (i = 0; i <= n-2; i++)
      f1.rep[i] = f.rep[n-i];
   f1.normalize();

   ZZ_pX f2;
   f2.rep.SetLength(n-1);
   for (i = 0; i <= n-2; i++)
      mul(f2.rep[i], f.rep[n-1-i], i+1);
   f2.normalize();

   ZZ_pX f3;
   InvTrunc(f3, f1, n-1);
   MulTrunc(f3, f3, f2, n-1);

   S.SetLength(n);

   S[0] = n;
   for (i = 1; i < n; i++)
      negate(S[i], coeff(f3, i-1));
}


void PlainTraceVec(vec_ZZ_p& S, const ZZ_pX& ff)
{
   if (deg(ff) <= 0)
      Error("TraceVec: bad args");

   ZZ_pX f;
   f = ff;

   MakeMonic(f);

   long n = deg(f);

   S.SetLength(n);

   if (n == 0)
      return;

   long k, i;
   ZZ acc, t;
   ZZ_p t1;

   S[0] = n;

   for (k = 1; k < n; k++) {
      mul(acc, rep(f.rep[n-k]), k);

      for (i = 1; i < k; i++) {
         mul(t, rep(f.rep[n-i]), rep(S[k-i]));
         add(acc, acc, t);
      }

      conv(t1, acc);
      negate(S[k], t1);
   }
}

void TraceVec(vec_ZZ_p& S, const ZZ_pX& f)
{
   if (deg(f) <= NTL_ZZ_pX_TRACE_CROSSOVER)
      PlainTraceVec(S, f);
   else
      FastTraceVec(S, f);
}

void ComputeTraceVec(const ZZ_pXModulus& F)
{
   vec_ZZ_p& S = *((vec_ZZ_p *) &F.tracevec);

   if (S.length() > 0)
      return;

   if (!F.UseFFT) {
      PlainTraceVec(S, F.f);
      return;
   }

   long i;
   long n = F.n;

   FFTRep R;
   ZZ_pX P, g;

   g.rep.SetLength(n-1);
   for (i = 1; i < n; i++)
      mul(g.rep[n-i-1], F.f.rep[n-i], i); 
   g.normalize();

   ToFFTRep(R, g, F.l);
   mul(R, R, F.HRep);
   FromFFTRep(P, R, n-2, 2*n-4);

   S.SetLength(n);

   S[0] = n;
   for (i = 1; i < n; i++)
      negate(S[i], coeff(P, n-1-i));
}

void TraceMod(ZZ_p& x, const ZZ_pX& a, const ZZ_pXModulus& F)
{
   long n = F.n;

   if (deg(a) >= n)
      Error("trace: bad args");

   if (F.tracevec.length() == 0) 
      ComputeTraceVec(F);

   InnerProduct(x, a.rep, F.tracevec);
}

void TraceMod(ZZ_p& x, const ZZ_pX& a, const ZZ_pX& f)
{
   if (deg(a) >= deg(f) || deg(f) <= 0)
      Error("trace: bad args");

   project(x, TraceVec(f), a);
}

void PlainResultant(ZZ_p& rres, const ZZ_pX& a, const ZZ_pX& b)
{
   ZZ_p res;
 
   if (IsZero(a) || IsZero(b))
      clear(res);
   else if (deg(a) == 0 && deg(b) == 0) 
      set(res);
   else {
      long d0, d1, d2;
      ZZ_p lc;
      set(res);

      long n = max(deg(a),deg(b)) + 1;
      ZZ_pX u(INIT_SIZE, n), v(INIT_SIZE, n);
      ZZVec tmp(n, ZZ_pInfo->ExtendedModulusSize);

      u = a;
      v = b;

      for (;;) {
         d0 = deg(u);
         d1 = deg(v);
         lc = LeadCoeff(v);

         PlainRem(u, u, v, tmp);
         swap(u, v);

         d2 = deg(v);
         if (d2 >= 0) {
            power(lc, lc, d0-d2);
            mul(res, res, lc);
            if (d0 & d1 & 1) negate(res, res);
         }
         else {
            if (d1 == 0) {
               power(lc, lc, d0);
               mul(res, res, lc);
            }
            else
               clear(res);
        
            break;
         }
      }

      rres = res;
   }
}


void ResIterHalfGCD(ZZ_pXMatrix& M_out, ZZ_pX& U, ZZ_pX& V, long d_red,
                    vec_ZZ_p& cvec, vec_long& dvec)
{
   M_out(0,0).SetMaxLength(d_red);
   M_out(0,1).SetMaxLength(d_red);
   M_out(1,0).SetMaxLength(d_red);
   M_out(1,1).SetMaxLength(d_red);

   set(M_out(0,0));   clear(M_out(0,1));
   clear(M_out(1,0)); set(M_out(1,1));

   long goal = deg(U) - d_red;

   if (deg(V) <= goal)
      return;

   ZZVec tmp(deg(U)+1, ZZ_pInfo->ExtendedModulusSize);
   ZZ_pX Q, t(INIT_SIZE, d_red);


   while (deg(V) > goal) {
      append(cvec, LeadCoeff(V));
      append(dvec, dvec[dvec.length()-1]-deg(U)+deg(V));
      PlainDivRem(Q, U, U, V, tmp);
      swap(U, V);

      mul(t, Q, M_out(1,0));
      sub(t, M_out(0,0), t);
      M_out(0,0) = M_out(1,0);
      M_out(1,0) = t;

      mul(t, Q, M_out(1,1));
      sub(t, M_out(0,1), t);
      M_out(0,1) = M_out(1,1);
      M_out(1,1) = t;
   }
}
   


void ResHalfGCD(ZZ_pXMatrix& M_out, const ZZ_pX& U, const ZZ_pX& V, long d_red,
                vec_ZZ_p& cvec, vec_long& dvec)
{
   if (IsZero(V) || deg(V) <= deg(U) - d_red) {
      set(M_out(0,0));   clear(M_out(0,1));
      clear(M_out(1,0)); set(M_out(1,1));
 
      return;
   }


   long n = deg(U) - 2*d_red + 2;
   if (n < 0) n = 0;

   ZZ_pX U1, V1;

   RightShift(U1, U, n);
   RightShift(V1, V, n);

   if (d_red <= NTL_ZZ_pX_HalfGCD_CROSSOVER) { 
      ResIterHalfGCD(M_out, U1, V1, d_red, cvec, dvec);
      return;
   }

   long d1 = (d_red + 1)/2;
   if (d1 < 1) d1 = 1;
   if (d1 >= d_red) d1 = d_red - 1;

   ZZ_pXMatrix M1;

   ResHalfGCD(M1, U1, V1, d1, cvec, dvec);
   mul(U1, V1, M1);

   long d2 = deg(V1) - deg(U) + n + d_red;

   if (IsZero(V1) || d2 <= 0) {
      M_out = M1;
      return;
   }


   ZZ_pX Q;
   ZZ_pXMatrix M2;

   append(cvec, LeadCoeff(V1));
   append(dvec, dvec[dvec.length()-1]-deg(U1)+deg(V1));
   DivRem(Q, U1, U1, V1);
   swap(U1, V1);

   ResHalfGCD(M2, U1, V1, d2, cvec, dvec);

   ZZ_pX t(INIT_SIZE, deg(M1(1,1))+deg(Q)+1);

   mul(t, Q, M1(1,0));
   sub(t, M1(0,0), t);
   swap(M1(0,0), M1(1,0));
   swap(M1(1,0), t);

   t.kill();

   t.SetMaxLength(deg(M1(1,1))+deg(Q)+1);

   mul(t, Q, M1(1,1));
   sub(t, M1(0,1), t);
   swap(M1(0,1), M1(1,1));
   swap(M1(1,1), t);

   t.kill();

   mul(M_out, M2, M1); 
}

void ResHalfGCD(ZZ_pX& U, ZZ_pX& V, vec_ZZ_p& cvec, vec_long& dvec)
{
   long d_red = (deg(U)+1)/2;

   if (IsZero(V) || deg(V) <= deg(U) - d_red) {
      return;
   }

   long du = deg(U);


   long d1 = (d_red + 1)/2;
   if (d1 < 1) d1 = 1;
   if (d1 >= d_red) d1 = d_red - 1;

   ZZ_pXMatrix M1;

   ResHalfGCD(M1, U, V, d1, cvec, dvec);
   mul(U, V, M1);

   long d2 = deg(V) - du + d_red;

   if (IsZero(V) || d2 <= 0) {
      return;
   }

   M1(0,0).kill();
   M1(0,1).kill();
   M1(1,0).kill();
   M1(1,1).kill();


   ZZ_pX Q;

   append(cvec, LeadCoeff(V));
   append(dvec, dvec[dvec.length()-1]-deg(U)+deg(V));
   DivRem(Q, U, U, V);
   swap(U, V);

   ResHalfGCD(M1, U, V, d2, cvec, dvec);

   mul(U, V, M1); 
}


void resultant(ZZ_p& rres, const ZZ_pX& u, const ZZ_pX& v)
{
   if (deg(u) <= NTL_ZZ_pX_GCD_CROSSOVER || deg(v) <= NTL_ZZ_pX_GCD_CROSSOVER) { 
      PlainResultant(rres, u, v);
      return;
   }

   ZZ_pX u1, v1;

   u1 = u;
   v1 = v;

   ZZ_p res, t;
   set(res);

   if (deg(u1) == deg(v1)) {
      rem(u1, u1, v1);
      swap(u1, v1);

      if (IsZero(v1)) {
         clear(rres);
         return;
      }

      power(t, LeadCoeff(u1), deg(u1) - deg(v1));
      mul(res, res, t);
      if (deg(u1) & 1)
         negate(res, res);
   }
   else if (deg(u1) < deg(v1)) {
      swap(u1, v1);
      if (deg(u1) & deg(v1) & 1)
         negate(res, res);
   }

   // deg(u1) > deg(v1) && v1 != 0

   vec_ZZ_p cvec;
   vec_long  dvec;

   cvec.SetMaxLength(deg(v1)+2);
   dvec.SetMaxLength(deg(v1)+2);

   append(cvec, LeadCoeff(u1));
   append(dvec, deg(u1));


   while (deg(u1) > NTL_ZZ_pX_GCD_CROSSOVER && !IsZero(v1)) { 
      ResHalfGCD(u1, v1, cvec, dvec);

      if (!IsZero(v1)) {
         append(cvec, LeadCoeff(v1));
         append(dvec, deg(v1));
         rem(u1, u1, v1);
         swap(u1, v1);
      }
   }

   if (IsZero(v1) && deg(u1) > 0) {
      clear(rres);
      return;
   }

   long i, l;
   l = dvec.length();

   if (deg(u1) == 0) {
      // we went all the way...

      for (i = 0; i <= l-3; i++) {
         power(t, cvec[i+1], dvec[i]-dvec[i+2]);
         mul(res, res, t);
         if (dvec[i] & dvec[i+1] & 1)
            negate(res, res);
      }

      power(t, cvec[l-1], dvec[l-2]);
      mul(res, res, t);
   }
   else {
      for (i = 0; i <= l-3; i++) {
         power(t, cvec[i+1], dvec[i]-dvec[i+2]);
         mul(res, res, t);
         if (dvec[i] & dvec[i+1] & 1)
            negate(res, res);
      }

      power(t, cvec[l-1], dvec[l-2]-deg(v1));
      mul(res, res, t);
      if (dvec[l-2] & dvec[l-1] & 1)
         negate(res, res);

      PlainResultant(t, u1, v1);
      mul(res, res, t);
   }

   rres = res;
}

void NormMod(ZZ_p& x, const ZZ_pX& a, const ZZ_pX& f)
{
   if (deg(f) <= 0 || deg(a) >= deg(f)) 
      Error("norm: bad args");

   if (IsZero(a)) {
      clear(x);
      return;
   }

   ZZ_p t;
   resultant(t, f, a);
   if (!IsOne(LeadCoeff(f))) {
      ZZ_p t1;
      power(t1, LeadCoeff(f), deg(a));
      inv(t1, t1);
      mul(t, t, t1);
   }

   x = t;
}

NTL_END_IMPL