lzz_pex.c

密码大家Shoup写的数论算法c语言实现

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

void IrredPolyMod(zz_pEX& h, const zz_pEX& g, const zz_pEXModulus& F, long m)
{
   if (m < 1 || m > F.n) Error("IrredPoly: bad args");

   zz_pEX R;
   set(R);

   DoMinPolyMod(h, g, F, m, R);
}



void IrredPolyMod(zz_pEX& h, const zz_pEX& g, const zz_pEXModulus& F)
{
   IrredPolyMod(h, g, F, F.n);
}



void MinPolyMod(zz_pEX& hh, const zz_pEX& g, const zz_pEXModulus& F)
{
   MinPolyMod(hh, g, F, F.n);
}

void diff(zz_pEX& x, const zz_pEX& 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_pEX& x)
{
   if (IsZero(x))
      return;

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

   zz_pE t;

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


long divide(zz_pEX& q, const zz_pEX& a, const zz_pEX& b)
{
   if (IsZero(b)) {
      if (IsZero(a)) {
         clear(q);
         return 1;
      }
      else
         return 0;
   }

   zz_pEX lq, r;
   DivRem(lq, r, a, b);
   if (!IsZero(r)) return 0; 
   q = lq;
   return 1;
}

long divide(const zz_pEX& a, const zz_pEX& b)
{
   if (IsZero(b)) return IsZero(a);
   zz_pEX lq, r;
   DivRem(lq, r, a, b);
   if (!IsZero(r)) return 0; 
   return 1;
}



static
long OptWinSize(long n)
// finds k that minimizes n/(k+1) + 2^{k-1}

{
   long k;
   double v, v_new;


   v = n/2.0 + 1.0;
   k = 1;

   for (;;) {
      v_new = n/(double(k+2)) + double(1L << k);
      if (v_new >= v) break;
      v = v_new;
      k++;
   }

   return k;
}
      


void PowerMod(zz_pEX& h, const zz_pEX& g, const ZZ& e, const zz_pEXModulus& F)
// h = g^e mod f using "sliding window" algorithm
{
   if (deg(g) >= F.n) Error("PowerMod: bad args");

   if (e == 0) {
      set(h);
      return;
   }

   if (e == 1) {
      h = g;
      return;
   }

   if (e == -1) {
      InvMod(h, g, F);
      return;
   }

   if (e == 2) {
      SqrMod(h, g, F);
      return;
   }

   if (e == -2) {
      SqrMod(h, g, F);
      InvMod(h, h, F);
      return;
   }


   long n = NumBits(e);

   zz_pEX res;
   res.SetMaxLength(F.n);
   set(res);

   long i;

   if (n < 16) {
      // plain square-and-multiply algorithm

      for (i = n - 1; i >= 0; i--) {
         SqrMod(res, res, F);
         if (bit(e, i))
            MulMod(res, res, g, F);
      }

      if (e < 0) InvMod(res, res, F);

      h = res;
      return;
   }

   long k = OptWinSize(n);
   k = min(k, 3);

   vec_zz_pEX v;

   v.SetLength(1L << (k-1));

   v[0] = g;
 
   if (k > 1) {
      zz_pEX t;
      SqrMod(t, g, F);

      for (i = 1; i < (1L << (k-1)); i++)
         MulMod(v[i], v[i-1], t, F);
   }


   long val;
   long cnt;
   long m;

   val = 0;
   for (i = n-1; i >= 0; i--) {
      val = (val << 1) | bit(e, i); 
      if (val == 0)
         SqrMod(res, res, F);
      else if (val >= (1L << (k-1)) || i == 0) {
         cnt = 0;
         while ((val & 1) == 0) {
            val = val >> 1;
            cnt++;
         }

         m = val;
         while (m > 0) {
            SqrMod(res, res, F);
            m = m >> 1;
         }

         MulMod(res, res, v[val >> 1], F);

         while (cnt > 0) {
            SqrMod(res, res, F);
            cnt--;
         }

         val = 0;
      }
   }

   if (e < 0) InvMod(res, res, F);

   h = res;
}

void InvMod(zz_pEX& x, const zz_pEX& a, const zz_pEX& f)
{
   if (deg(a) >= deg(f) || deg(f) == 0) Error("InvMod: bad args");

   zz_pEX d, t;

   XGCD(d, x, t, a, f);
   if (!IsOne(d))
      Error("zz_pEX InvMod: can't compute multiplicative inverse");
}

long InvModStatus(zz_pEX& x, const zz_pEX& a, const zz_pEX& f)
{
   if (deg(a) >= deg(f) || deg(f) == 0) Error("InvModStatus: bad args");
   zz_pEX d, t;

   XGCD(d, x, t, a, f);
   if (!IsOne(d)) {
      x = d;
      return 1;
   }
   else
      return 0;
}


void MulMod(zz_pEX& x, const zz_pEX& a, const zz_pEX& b, const zz_pEX& f)
{
   if (deg(a) >= deg(f) || deg(b) >= deg(f) || deg(f) == 0)
      Error("MulMod: bad args");

   zz_pEX t;

   mul(t, a, b);
   rem(x, t, f);
}

void SqrMod(zz_pEX& x, const zz_pEX& a, const zz_pEX& f)
{
   if (deg(a) >= deg(f) || deg(f) == 0) Error("SqrMod: bad args");

   zz_pEX t;

   sqr(t, a);
   rem(x, t, f);
}


void PowerXMod(zz_pEX& hh, const ZZ& e, const zz_pEXModulus& F)
{
   if (F.n < 0) Error("PowerXMod: uninitialized modulus");

   if (IsZero(e)) {
      set(hh);
      return;
   }

   long n = NumBits(e);
   long i;

   zz_pEX h;

   h.SetMaxLength(F.n);
   set(h);

   for (i = n - 1; i >= 0; i--) {
      SqrMod(h, h, F);
      if (bit(e, i))
         MulByXMod(h, h, F.f);
   }

   if (e < 0) InvMod(h, h, F);

   hh = h;
}


void reverse(zz_pEX& x, const zz_pEX& a, long hi)
{
   if (hi < -1) Error("reverse: bad args");
   if (hi >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("overflow in reverse");

   if (&x == &a) {
      zz_pEX tmp;
      CopyReverse(tmp, a, hi);
      x = tmp;
   }
   else
      CopyReverse(x, a, hi);
}


void power(zz_pEX& x, const zz_pEX& a, long e)
{
   if (e < 0) {
      Error("power: negative exponent");
   }

   if (e == 0) {
      x = 1;
      return;
   }

   if (a == 0 || a == 1) {
      x = a;
      return;
   }

   long da = deg(a);

   if (da == 0) {
      x = power(ConstTerm(a), e);
      return;
   }

   if (da > (NTL_MAX_LONG-1)/e)
      Error("overflow in power");

   zz_pEX res;
   res.SetMaxLength(da*e + 1);
   res = 1;
   
   long k = NumBits(e);
   long i;

   for (i = k - 1; i >= 0; i--) {
      sqr(res, res);
      if (bit(e, i))
         mul(res, res, a);
   }

   x = res;
}



static
void FastTraceVec(vec_zz_pE& S, const zz_pEXModulus& f)
{
   long n = deg(f);

   zz_pEX x = reverse(-LeftShift(reverse(diff(reverse(f)), n-1), n-1)/f, n-1);

   S.SetLength(n);
   S[0] = n;

   long i;
   for (i = 1; i < n; i++)
      S[i] = coeff(x, i);
}


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

   zz_pEX f;
   f = ff;

   MakeMonic(f);

   long n = deg(f);

   S.SetLength(n);

   if (n == 0)
      return;

   long k, i;
   zz_pX acc, t;
   zz_pE 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_pE& S, const zz_pEX& f)
{
   if (deg(f) < zz_pE::DivCross())
      PlainTraceVec(S, f);
   else
      FastTraceVec(S, f);
}

static
void ComputeTraceVec(const zz_pEXModulus& F)
{
   vec_zz_pE& S = *((vec_zz_pE *) &F.tracevec);

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

   if (F.method == zz_pEX_MOD_PLAIN) {
      PlainTraceVec(S, F.f);
   }
   else {
      FastTraceVec(S, F);
   }
}

void TraceMod(zz_pE& x, const zz_pEX& a, const zz_pEXModulus& 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_pE& x, const zz_pEX& a, const zz_pEX& f)
{
   if (deg(a) >= deg(f) || deg(f) <= 0)
      Error("trace: bad args");

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


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

      long n = max(deg(a),deg(b)) + 1;
      zz_pEX u(INIT_SIZE, n), v(INIT_SIZE, n);
      vec_zz_pX tmp;
      SetSize(tmp, n, 2*zz_pE::degree());

      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 resultant(zz_pE& rres, const zz_pEX& a, const zz_pEX& b)
{
   PlainResultant(rres, a, b); 
}


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

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

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

   x = t;
}



// tower stuff...



void InnerProduct(zz_pEX& x, const vec_zz_p& v, long low, long high,
                   const vec_zz_pEX& H, long n, vec_zz_pE& t)
{
   zz_pE s;
   long i, j;

   for (j = 0; j < n; j++)
      clear(t[j]);

   high = min(high, v.length()-1);
   for (i = low; i <= high; i++) {
      const vec_zz_pE& h = H[i-low].rep;
      long m = h.length();
      const zz_p& w = v[i];

      for (j = 0; j < m; j++) {
         mul(s, h[j], w);
         add(t[j], t[j], s);
      }
   }

   x.rep.SetLength(n);
   for (j = 0; j < n; j++)
      x.rep[j] = t[j];

   x.normalize();
}



void CompTower(zz_pEX& x, const zz_pX& g, const zz_pEXArgument& A,
             const zz_pEXModulus& F)
{
   if (deg(g) <= 0) {
      conv(x, g);
      return;
   }


   zz_pEX s, t;
   vec_zz_pE scratch;
   scratch.SetLength(deg(F));

   long m = A.H.length() - 1;
   long l = ((g.rep.length()+m-1)/m) - 1;

   const zz_pEX& M = A.H[m];

   InnerProduct(t, g.rep, l*m, l*m + m - 1, A.H, F.n, scratch);
   for (long i = l-1; i >= 0; i--) {
      InnerProduct(s, g.rep, i*m, i*m + m - 1, A.H, F.n, scratch);
      MulMod(t, t, M, F);
      add(t, t, s);
   }
   x = t;
}


void CompTower(zz_pEX& x, const zz_pX& g, const zz_pEX& h, 
             const zz_pEXModulus& F)
   // x = g(h) mod f
{
   long m = SqrRoot(g.rep.length());

   if (m == 0) {
      clear(x);
      return;
   }


   zz_pEXArgument A;

   build(A, h, F, m);

   CompTower(x, g, A, F);
}

void PrepareProjection(vec_vec_zz_p& tt, const vec_zz_pE& s,
                       const vec_zz_p& proj)
{
   long l = s.length();
   tt.SetLength(l);

   zz_pXMultiplier M;
   long i;

   for (i = 0; i < l; i++) {
      build(M, rep(s[i]), zz_pE::modulus());
      UpdateMap(tt[i], proj, M, zz_pE::modulus());
   }
}

void ProjectedInnerProduct(zz_p& x, const vec_zz_pE& a, 
                           const vec_vec_zz_p& b)
{
   long n = min(a.length(), b.length());

   zz_p t, res;

   res = 0;

   long i;
   for (i = 0; i < n; i++) {
      project(t, b[i], rep(a[i]));
      res += t;
   }

   x = res;
}


   
void PrecomputeProj(vec_zz_p& proj, const zz_pX& f)
{
   long n = deg(f);

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

   if (ConstTerm(f) != 0) {
      proj.SetLength(1);
      proj[0] = 1;
   }
   else {
      proj.SetLength(n);
      clear(proj);
      proj[n-1] = 1;
   }
}


void ProjectPowersTower(vec_zz_p& x, const vec_zz_pE& a, long k,
                   const zz_pEXArgument& H, const zz_pEXModulus& F,
                   const vec_zz_p& proj)

{
   long n = F.n;

   if (a.length() > n || k < 0 || k >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("ProjectPowers: bad args");

   long m = H.H.length()-1;
   long l = (k+m-1)/m - 1;

   zz_pEXTransMultiplier M;
   build(M, H.H[m], F);

   vec_zz_pE s(INIT_SIZE, n);
   s = a;

   x.SetLength(k);

   vec_vec_zz_p tt;

   for (long i = 0; i <= l; i++) {
      long m1 = min(m, k-i*m);
      zz_p* w = &x[i*m];

      PrepareProjection(tt, s, proj);

      for (long j = 0; j < m1; j++)
         ProjectedInnerProduct(w[j], H.H[j].rep, tt);
      if (i < l)
         UpdateMap(s, s, M, F);
   }
}




void ProjectPowersTower(vec_zz_p& x, const vec_zz_pE& a, long k,
                   const zz_pEX& h, const zz_pEXModulus& F,
                   const vec_zz_p& proj)

{
   if (a.length() > F.n || k < 0) Error("ProjectPowers: bad args");

   if (k == 0) {
      x.SetLength(0);
      return;
   }

   long m = SqrRoot(k);

   zz_pEXArgument H;

   build(H, h, F, m);
   ProjectPowersTower(x, a, k, H, F, proj);
}


void DoMinPolyTower(zz_pX& h, const zz_pEX& g, const zz_pEXModulus& F, long m,
               const vec_zz_pE& R, const vec_zz_p& proj)
{
   vec_zz_p x;

   ProjectPowersTower(x, R, 2*m, g, F, proj);
   
   MinPolySeq(h, x, m);
}


void ProbMinPolyTower(zz_pX& h, const zz_pEX& g, const zz_pEXModulus& F, 
                      long m)
{
   long n = F.n;
   if (m < 1 || m > n*zz_pE::degree()) Error("ProbMinPoly: bad args");

   vec_zz_pE R;
   R.SetLength(n);
   long i;
   for (i = 0; i < n; i++)
      random(R[i]);

   vec_zz_p proj;
   PrecomputeProj(proj, zz_pE::modulus());

   DoMinPolyTower(h, g, F, m, R, proj);
}


void ProbMinPolyTower(zz_pX& h, const zz_pEX& g, const zz_pEXModulus& F, 
                      long m, const vec_zz_p& proj)
{
   long n = F.n;
   if (m < 1 || m > n*zz_pE::degree()) Error("ProbMinPoly: bad args");

   vec_zz_pE R;
   R.SetLength(n);
   long i;
   for (i = 0; i < n; i++)
      random(R[i]);

   DoMinPolyTower(h, g, F, m, R, proj);
}

void MinPolyTower(zz_pX& hh, const zz_pEX& g, const zz_pEXModulus& F, long m)
{
   zz_pX h;
   zz_pEX h1;
   long n = F.n;
   if (m < 1 || m > n*zz_pE::degree()) {
      Error("MinPoly: bad args");
   }

   vec_zz_p proj;
   PrecomputeProj(proj, zz_pE::modulus());

   /* probabilistically compute min-poly */

   ProbMinPolyTower(h, g, F, m, proj);
   if (deg(h) == m) { hh = h; return; }
   CompTower(h1, h, g, F);
   if (IsZero(h1)) { hh = h; return; }

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

   long i;

   zz_pX h2;
   zz_pEX h3;
   vec_zz_pE R;
   zz_pEXTransMultiplier H1;
   

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

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

void IrredPolyTower(zz_pX& h, const zz_pEX& g, const zz_pEXModulus& F, long m)
{
   if (m < 1 || m > deg(F)*zz_pE::degree()) Error("IrredPoly: bad args");

   vec_zz_pE R;
   R.SetLength(1);
   R[0] = 1;

   vec_zz_p proj;
   proj.SetLength(1);
   proj[0] = 1;

   DoMinPolyTower(h, g, F, m, R, proj);
}

NTL_END_IMPL