lzz_pex.c

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

   if (deg(a) >= F.n || deg(b) >= F.n) Error("MulMod: bad args");

   zz_pEX t;
   mul(t, a, b);
   rem(c, t, F);
}


void SqrMod(zz_pEX& c, const zz_pEX& a, const zz_pEXModulus& F)
{
   if (deg(a) >= F.n) Error("MulMod: bad args");

   zz_pEX t;
   sqr(t, a);
   rem(c, t, F);
}



void UseMulRem(zz_pEX& r, const zz_pEX& a, const zz_pEX& b)
{
   zz_pEX P1;
   zz_pEX P2;

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

   CopyReverse(P1, b, db);
   InvTrunc(P2, P1, da-db+1);
   CopyReverse(P1, P2, da-db);

   RightShift(P2, a, db);
   mul(P2, P1, P2);
   RightShift(P2, P2, da-db);
   mul(P1, P2, b);
   sub(P1, a, P1);
   
   r = P1;
}

void UseMulDivRem(zz_pEX& q, zz_pEX& r, const zz_pEX& a, const zz_pEX& b)
{
   zz_pEX P1;
   zz_pEX P2;

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

   CopyReverse(P1, b, db);
   InvTrunc(P2, P1, da-db+1);
   CopyReverse(P1, P2, da-db);

   RightShift(P2, a, db);
   mul(P2, P1, P2);
   RightShift(P2, P2, da-db);
   mul(P1, P2, b);
   sub(P1, a, P1);
   
   r = P1;
   q = P2;
}

void UseMulDiv(zz_pEX& q, const zz_pEX& a, const zz_pEX& b)
{
   zz_pEX P1;
   zz_pEX P2;

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

   CopyReverse(P1, b, db);
   InvTrunc(P2, P1, da-db+1);
   CopyReverse(P1, P2, da-db);

   RightShift(P2, a, db);
   mul(P2, P1, P2);
   RightShift(P2, P2, da-db);
   
   q = P2;
}



void DivRem(zz_pEX& q, zz_pEX& r, const zz_pEX& a, const zz_pEX& b)
{
   long sa = a.rep.length();
   long sb = b.rep.length();

   if (sb < zz_pE::DivCross() || sa-sb < zz_pE::DivCross())
      PlainDivRem(q, r, a, b);
   else if (sa < 4*sb)
      UseMulDivRem(q, r, a, b);
   else {
      zz_pEXModulus B;
      build(B, b);
      DivRem(q, r, a, B);
   }
}

void div(zz_pEX& q, const zz_pEX& a, const zz_pEX& b)
{
   long sa = a.rep.length();
   long sb = b.rep.length();

   if (sb < zz_pE::DivCross() || sa-sb < zz_pE::DivCross())
      PlainDiv(q, a, b);
   else if (sa < 4*sb)
      UseMulDiv(q, a, b);
   else {
      zz_pEXModulus B;
      build(B, b);
      div(q, a, B);
   }
}

void div(zz_pEX& q, const zz_pEX& a, const zz_pE& b)
{
   zz_pE T;
   inv(T, b);
   mul(q, a, T);
}

void div(zz_pEX& q, const zz_pEX& a, const zz_p& b)
{
   NTL_zz_pRegister(T);
   inv(T, b);
   mul(q, a, T);
}

void div(zz_pEX& q, const zz_pEX& a, long b)
{
   NTL_zz_pRegister(T);
   T = b;
   inv(T, T);
   mul(q, a, T);
}

void rem(zz_pEX& r, const zz_pEX& a, const zz_pEX& b)
{
   long sa = a.rep.length();
   long sb = b.rep.length();

   if (sb < zz_pE::DivCross() || sa-sb < zz_pE::DivCross())
      PlainRem(r, a, b);
   else if (sa < 4*sb)
      UseMulRem(r, a, b);
   else {
      zz_pEXModulus B;
      build(B, b);
      rem(r, a, B);
   }
}

void GCD(zz_pEX& x, const zz_pEX& a, const zz_pEX& b)
{
   zz_pE t;

   if (IsZero(b))
      x = a;
   else if (IsZero(a))
      x = b;
   else {
      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;
      do {
         PlainRem(u, u, v, tmp);
         swap(u, v);
      } while (!IsZero(v));

      x = u;
   }

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

   /* make gcd monic */


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



         

void XGCD(zz_pEX& d, zz_pEX& s, zz_pEX& t, const zz_pEX& a, const zz_pEX& b)
{
   zz_pE z;


   if (IsZero(b)) {
      set(s);
      clear(t);
      d = a;
   }
   else if (IsZero(a)) {
      clear(s);
      set(t);
      d = b;
   }
   else {
      long e = max(deg(a), deg(b)) + 1;

      zz_pEX temp(INIT_SIZE, e), u(INIT_SIZE, e), v(INIT_SIZE, e), 
            u0(INIT_SIZE, e), v0(INIT_SIZE, e), 
            u1(INIT_SIZE, e), v1(INIT_SIZE, e), 
            u2(INIT_SIZE, e), v2(INIT_SIZE, e), q(INIT_SIZE, e);


      set(u1); clear(v1);
      clear(u2); set(v2);
      u = a; v = b;

      do {
         DivRem(q, u, u, v);
         swap(u, v);
         u0 = u2;
         v0 = v2;
         mul(temp, q, u2);
         sub(u2, u1, temp);
         mul(temp, q, v2);
         sub(v2, v1, temp);
         u1 = u0;
         v1 = v0;
      } while (!IsZero(v));

      d = u;
      s = u1;
      t = v1;
   }

   if (IsZero(d)) return;
   if (IsOne(LeadCoeff(d))) return;

   /* make gcd monic */

   inv(z, LeadCoeff(d));
   mul(d, d, z);
   mul(s, s, z);
   mul(t, t, z);
}

NTL_vector_impl(zz_pEX,vec_zz_pEX)

NTL_eq_vector_impl(zz_pEX,vec_zz_pEX)

NTL_io_vector_impl(zz_pEX,vec_zz_pEX)

void IterBuild(zz_pE* a, long n)
{
   long i, k;
   zz_pE b, t;

   if (n <= 0) return;

   negate(a[0], a[0]);

   for (k = 1; k <= n-1; k++) {
      negate(b, a[k]);
      add(a[k], b, a[k-1]);
      for (i = k-1; i >= 1; i--) {
         mul(t, a[i], b);
         add(a[i], t, a[i-1]);
      }
      mul(a[0], a[0], b);
   }
}

void BuildFromRoots(zz_pEX& x, const vec_zz_pE& a)
{
   long n = a.length();

   if (n == 0) {
      set(x);
      return;
   }

   x.rep.SetMaxLength(n+1);
   x.rep = a;
   IterBuild(&x.rep[0], n);
   x.rep.SetLength(n+1);
   SetCoeff(x, n);
}

void eval(zz_pE& b, const zz_pEX& f, const zz_pE& a)
// does a Horner evaluation
{
   zz_pE acc;
   long i;

   clear(acc);
   for (i = deg(f); i >= 0; i--) {
      mul(acc, acc, a);
      add(acc, acc, f.rep[i]);
   }

   b = acc;
}

void eval(vec_zz_pE& b, const zz_pEX& f, const vec_zz_pE& a)
// naive algorithm:  repeats Horner
{
   if (&b == &f.rep) {
      vec_zz_pE bb;
      eval(bb, f, a);
      b = bb;
      return;
   }

   long m = a.length();
   b.SetLength(m);
   long i;
   for (i = 0; i < m; i++)
      eval(b[i], f, a[i]);
}


void interpolate(zz_pEX& f, const vec_zz_pE& a, const vec_zz_pE& b)
{
   long m = a.length();
   if (b.length() != m) Error("interpolate: vector length mismatch");

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

   vec_zz_pE prod;
   prod = a;

   zz_pE t1, t2;

   long k, i;

   vec_zz_pE res;
   res.SetLength(m);

   for (k = 0; k < m; k++) {

      const zz_pE& aa = a[k];

      set(t1);
      for (i = k-1; i >= 0; i--) {
         mul(t1, t1, aa);
         add(t1, t1, prod[i]);
      }

      clear(t2);
      for (i = k-1; i >= 0; i--) {
         mul(t2, t2, aa);
         add(t2, t2, res[i]);
      }


      inv(t1, t1);
      sub(t2, b[k], t2);
      mul(t1, t1, t2);

      for (i = 0; i < k; i++) {
         mul(t2, prod[i], t1);
         add(res[i], res[i], t2);
      }

      res[k] = t1;

      if (k < m-1) {
         if (k == 0)
            negate(prod[0], prod[0]);
         else {
            negate(t1, a[k]);
            add(prod[k], t1, prod[k-1]);
            for (i = k-1; i >= 1; i--) {
               mul(t2, prod[i], t1);
               add(prod[i], t2, prod[i-1]);
            }
            mul(prod[0], prod[0], t1);
         }
      }
   }

   while (m > 0 && IsZero(res[m-1])) m--;
   res.SetLength(m);
   f.rep = res;
}
   
void InnerProduct(zz_pEX& x, const vec_zz_pE& v, long low, long high, 
                   const vec_zz_pEX& H, long n, vec_zz_pX& t)
{
   zz_pX 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_pX& w = rep(v[i]);

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

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



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


   zz_pEX s, t;
   vec_zz_pX scratch;
   SetSize(scratch, deg(F), 2*zz_pE::degree());

   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 build(zz_pEXArgument& A, const zz_pEX& h, const zz_pEXModulus& F, long m)
{
   long i;

   if (m <= 0 || deg(h) >= F.n)
      Error("build: bad args");

   if (m > F.n) m = F.n;

   if (zz_pEXArgBound > 0) {
      double sz = zz_p::storage();
      sz = sz*zz_pE::degree();
      sz = sz + NTL_VECTOR_HEADER_SIZE + sizeof(vec_zz_p);
      sz = sz*F.n;
      sz = sz + NTL_VECTOR_HEADER_SIZE + sizeof(vec_zz_pE);
      sz = sz/1024;
      m = min(m, long(zz_pEXArgBound/sz));
      m = max(m, 1);
   }



   A.H.SetLength(m+1);

   set(A.H[0]);
   A.H[1] = h;
   for (i = 2; i <= m; i++)
      MulMod(A.H[i], A.H[i-1], h, F);
}

long zz_pEXArgBound = 0;




void CompMod(zz_pEX& x, const zz_pEX& 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);

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




void Comp2Mod(zz_pEX& x1, zz_pEX& x2, const zz_pEX& g1, const zz_pEX& g2,
              const zz_pEX& h, const zz_pEXModulus& F)

{
   long m = SqrRoot(g1.rep.length() + g2.rep.length());

   if (m == 0) {
      clear(x1);
      clear(x2);
      return;
   }

   zz_pEXArgument A;

   build(A, h, F, m);

   zz_pEX xx1, xx2;

   CompMod(xx1, g1, A, F);
   CompMod(xx2, g2, A, F);

   x1 = xx1;
   x2 = xx2;
}

void Comp3Mod(zz_pEX& x1, zz_pEX& x2, zz_pEX& x3, 
              const zz_pEX& g1, const zz_pEX& g2, const zz_pEX& g3,
              const zz_pEX& h, const zz_pEXModulus& F)

{
   long m = SqrRoot(g1.rep.length() + g2.rep.length() + g3.rep.length());

   if (m == 0) {
      clear(x1);
      clear(x2);
      clear(x3);
      return;
   }

   zz_pEXArgument A;

   build(A, h, F, m);

   zz_pEX xx1, xx2, xx3;

   CompMod(xx1, g1, A, F);
   CompMod(xx2, g2, A, F);
   CompMod(xx3, g3, A, F);

   x1 = xx1;
   x2 = xx2;
   x3 = xx3;
}

void build(zz_pEXTransMultiplier& B, const zz_pEX& b, const zz_pEXModulus& F)
{
   long db = deg(b);

   if (db >= F.n) Error("build TransMultiplier: bad args");

   zz_pEX t;

   LeftShift(t, b, F.n-1);
   div(t, t, F);

   // we optimize for low degree b

   long d;

   d = deg(t);
   if (d < 0)
      B.shamt_fbi = 0;
   else
      B.shamt_fbi = F.n-2 - d; 

   CopyReverse(B.fbi, t, d);

   // The following code optimizes the case when 
   // f = X^n + low degree poly

   trunc(t, F.f, F.n);
   d = deg(t);
   if (d < 0)
      B.shamt = 0;
   else
      B.shamt = d;

   CopyReverse(B.f0, t, d);

   if (db < 0)
      B.shamt_b = 0;
   else
      B.shamt_b = db;

   CopyReverse(B.b, b, db);
}

void TransMulMod(zz_pEX& x, const zz_pEX& a, const zz_pEXTransMultiplier& B,
               const zz_pEXModulus& F)
{
   if (deg(a) >= F.n) Error("TransMulMod: bad args");

   zz_pEX t1, t2;

   mul(t1, a, B.b);
   RightShift(t1, t1, B.shamt_b);

   mul(t2, a, B.f0);
   RightShift(t2, t2, B.shamt);
   trunc(t2, t2, F.n-1);

   mul(t2, t2, B.fbi);
   if (B.shamt_fbi > 0) LeftShift(t2, t2, B.shamt_fbi);
   trunc(t2, t2, F.n-1);
   LeftShift(t2, t2, 1);

   sub(x, t1, t2);
}


void ShiftSub(zz_pEX& U, const zz_pEX& V, long n)
// assumes input does not alias output
{
   if (IsZero(V))
      return;

   long du = deg(U);
   long dv = deg(V);

   long d = max(du, n+dv);

   U.rep.SetLength(d+1);
   long i;

   for (i = du+1; i <= d; i++)
      clear(U.rep[i]);

   for (i = 0; i <= dv; i++)
      sub(U.rep[i+n], U.rep[i+n], V.rep[i]);

   U.normalize();
}


void UpdateMap(vec_zz_pE& x, const vec_zz_pE& a,
         const zz_pEXTransMultiplier& B, const zz_pEXModulus& F)
{
   zz_pEX xx;
   TransMulMod(xx, to_zz_pEX(a), B, F);
   x = xx.rep;
}

static
void ProjectPowers(vec_zz_pE& x, const zz_pEX& a, long k, 
                   const zz_pEXArgument& H, const zz_pEXModulus& F)
{
   if (k < 0 || k >= (1L << (NTL_BITS_PER_LONG-4)) || deg(a) >= F.n)
      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);

   zz_pEX s;
   s = a;

   x.SetLength(k);

   long i;

   for (i = 0; i <= l; i++) {
      long m1 = min(m, k-i*m);
      for (long j = 0; j < m1; j++)
         InnerProduct(x[i*m+j], H.H[j].rep, s.rep);
      if (i < l)
         TransMulMod(s, s, M, F);
   }
}

static
void ProjectPowers(vec_zz_pE& x, const zz_pEX& a, long k, const zz_pEX& h, 
                   const zz_pEXModulus& F)
{
   if (k < 0 || deg(a) >= F.n || deg(h) >= F.n)
      Error("ProjectPowers: bad args");

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

   long m = SqrRoot(k);

   zz_pEXArgument H;
   build(H, h, F, m);

   ProjectPowers(x, a, k, H, F);
}

void ProjectPowers(vec_zz_pE& x, const vec_zz_pE& a, long k,
                   const zz_pEXArgument& H, const zz_pEXModulus& F)
{
   ProjectPowers(x, to_zz_pEX(a), k, H, F);
}

void ProjectPowers(vec_zz_pE& x, const vec_zz_pE& a, long k,
                   const zz_pEX& h, const zz_pEXModulus& F)
{
   ProjectPowers(x, to_zz_pEX(a), k, h, F);
}




void BerlekampMassey(zz_pEX& h, const vec_zz_pE& a, long m)
{
   zz_pEX Lambda, Sigma, Temp;
   long L;
   zz_pE Delta, Delta1, t1;
   long shamt;

   // cerr << "*** " << m << "\n";

   Lambda.SetMaxLength(m+1);
   Sigma.SetMaxLength(m+1);
   Temp.SetMaxLength(m+1);

   L = 0;
   set(Lambda);
   clear(Sigma);
   set(Delta);
   shamt = 0;

   long i, r, dl;

   for (r = 1; r <= 2*m; r++) {
      // cerr << r << "--";
      clear(Delta1);
      dl = deg(Lambda);
      for (i = 0; i <= dl; i++) {
         mul(t1, Lambda.rep[i], a[r-i-1]);
         add(Delta1, Delta1, t1);
      }

      if (IsZero(Delta1)) {
         shamt++;
         // cerr << "case 1: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
      }
      else if (2*L < r) {
         div(t1, Delta1, Delta);
         mul(Temp, Sigma, t1);
         Sigma = Lambda;
         ShiftSub(Lambda, Temp, shamt+1);
         shamt = 0;
         L = r-L;
         Delta = Delta1;
         // cerr << "case 2: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
      }
      else {
         shamt++;
         div(t1, Delta1, Delta);
         mul(Temp, Sigma, t1);
         ShiftSub(Lambda, Temp, shamt);
         // cerr << "case 3: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
      }
   }

   // cerr << "finished: " << L << " " << deg(Lambda) << "\n"; 

   dl = deg(Lambda);
   h.rep.SetLength(L + 1);

   for (i = 0; i < L - dl; i++)
      clear(h.rep[i]);

   for (i = L - dl; i <= L; i++)
      h.rep[i] = Lambda.rep[L - i];
}




void MinPolySeq(zz_pEX& h, const vec_zz_pE& a, long m)
{
   if (m < 0 || m >= (1L << (NTL_BITS_PER_LONG-4))) Error("MinPoly: bad args");
   if (a.length() < 2*m) Error("BerlekampMassey: sequence too short");

   BerlekampMassey(h, a, m);
}


void DoMinPolyMod(zz_pEX& h, const zz_pEX& g, const zz_pEXModulus& F, long m, 
               const zz_pEX& R)
{
   vec_zz_pE x;

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

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

   zz_pEX R;
   random(R, n);

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

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

void MinPolyMod(zz_pEX& hh, const zz_pEX& g, const zz_pEXModulus& F, long m)
{
   zz_pEX 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 */

   zz_pEX h2, h3;
   zz_pEX R;
   zz_pEXTransMultiplier H1;
   

   for (;;) {
      random(R, n);
      build(H1, h1, F);
      TransMulMod(R, R, H1, F);
      DoMinPolyMod(h2, g, F, m-deg(h), R);

      mul(h, h, h2);
      if (deg(h) == m) { hh = h; return; }