zzx1.c

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

   ap = a.rep.elts();
   xp = x.rep.elts();

   for (i = 0; i <= da; i++) 
      mul(xp[i], ap[i], b);
}




void diff(ZZX& x, const ZZX& 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 HomPseudoDivRem(ZZX& q, ZZX& r, const ZZX& a, const ZZX& b)
{
   if (IsZero(b)) Error("division by zero");

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

   if (da < db) {
      r = b;
      clear(q);
      return;
   }

   ZZ LC;
   LC = LeadCoeff(b);

   ZZ LC1;

   power(LC1, LC, da-db+1);

   long a_bound = NumBits(LC1) + MaxBits(a);

   LC1.kill();

   long b_bound = MaxBits(b);

   zz_pBak bak;
   bak.save();

   ZZX qq, rr;

   ZZ prod, t;
   set(prod);

   clear(qq);
   clear(rr);

   long i;
   long Qinstable, Rinstable;

   Qinstable = 1;
   Rinstable = 1;

   for (i = 0; ; i++) {
      zz_p::FFTInit(i);
      long p = zz_p::modulus();


      if (divide(LC, p)) continue;

      zz_pX A, B, Q, R;

      conv(A, a);
      conv(B, b);
      
      if (!IsOne(LC)) {
         zz_p y;
         conv(y, LC);
         power(y, y, da-db+1);
         mul(A, A, y);
      }

      if (!Qinstable) {
         conv(Q, qq);
         mul(R, B, Q);
         sub(R, A, R);

         if (deg(R) >= db)
            Qinstable = 1;
         else
            Rinstable = CRT(rr, prod, R);
      }

      if (Qinstable) {
         DivRem(Q, R, A, B);
         t = prod;
         Qinstable = CRT(qq, t, Q);
         Rinstable =  CRT(rr, prod, R);
      }

      if (!Qinstable && !Rinstable) {
         // stabilized...check if prod is big enough

         long bound1 = b_bound + MaxBits(qq) + NumBits(min(db, da-db)+1);
         long bound2 = MaxBits(rr);
         long bound = max(bound1, bound2);

         if (a_bound > bound)
            bound = a_bound;

         bound += 4;

         if (NumBits(prod) > bound)
            break;
      }
   }

   bak.restore();

   q = qq;
   r = rr;
}




void HomPseudoDiv(ZZX& q, const ZZX& a, const ZZX& b)
{
   ZZX r;
   HomPseudoDivRem(q, r, a, b);
}

void HomPseudoRem(ZZX& r, const ZZX& a, const ZZX& b)
{
   ZZX q;
   HomPseudoDivRem(q, r, a, b);
}

void PlainPseudoDivRem(ZZX& q, ZZX& r, const ZZX& a, const ZZX& b)
{
   long da, db, dq, i, j, LCIsOne;
   const ZZ *bp;
   ZZ *qp;
   ZZ *xp;


   ZZ  s, t;

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

   if (db < 0) Error("ZZX: division by zero");

   if (da < db) {
      r = a;
      clear(q);
      return;
   }

   ZZX lb;

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

   ZZ LC = bp[db];
   LCIsOne = IsOne(LC);


   vec_ZZ x;

   x = a.rep;
   xp = x.elts();

   dq = da - db;
   q.rep.SetLength(dq+1);
   qp = q.rep.elts();

   if (!LCIsOne) {
      t = LC;
      for (i = dq-1; i >= 0; i--) {
         mul(xp[i], xp[i], t);
         if (i > 0) mul(t, t, LC);
      }
   }

   for (i = dq; i >= 0; i--) {
      t = xp[i+db];
      qp[i] = t;

      for (j = db-1; j >= 0; j--) {
	 mul(s, t, bp[j]);
         if (!LCIsOne) mul(xp[i+j], xp[i+j], LC);
	 sub(xp[i+j], xp[i+j], s);
      }
   }

   if (!LCIsOne) {
      t = LC;
      for (i = 1; i <= dq; i++) {
         mul(qp[i], qp[i], t);
         if (i < dq) mul(t, t, LC);
      }
   }
      

   r.rep.SetLength(db);
   for (i = 0; i < db; i++)
      r.rep[i] = xp[i];
   r.normalize();
}


void PlainPseudoDiv(ZZX& q, const ZZX& a, const ZZX& b)
{
   ZZX r;
   PlainPseudoDivRem(q, r, a, b);
}

void PlainPseudoRem(ZZX& r, const ZZX& a, const ZZX& b)
{
   ZZX q;
   PlainPseudoDivRem(q, r, a, b);
}

void div(ZZX& q, const ZZX& a, long b)
{
   if (b == 0) Error("div: division by zero");

   if (!divide(q, a, b)) Error("DivRem: quotient undefined over ZZ");
}

void div(ZZX& q, const ZZX& a, const ZZ& b)
{
   if (b == 0) Error("div: division by zero");

   if (!divide(q, a, b)) Error("DivRem: quotient undefined over ZZ");
}

static
void ConstDivRem(ZZX& q, ZZX& r, const ZZX& a, const ZZ& b)
{
   if (b == 0) Error("DivRem: division by zero");

   if (!divide(q, a, b)) Error("DivRem: quotient undefined over ZZ");

   r = 0;
}

static
void ConstRem(ZZX& r, const ZZX& a, const ZZ& b)
{
   if (b == 0) Error("rem: division by zero");

   r = 0;
}



void DivRem(ZZX& q, ZZX& r, const ZZX& a, const ZZX& b)
{
   long da = deg(a);
   long db = deg(b);

   if (db < 0) Error("DivRem: division by zero");

   if (da < db) {
      r = a;
      q = 0;
   }
   else if (db == 0) {
      ConstDivRem(q, r, a, ConstTerm(b));
   }
   else if (IsOne(LeadCoeff(b))) {
      PseudoDivRem(q, r, a, b);
   }
   else if (LeadCoeff(b) == -1) {
      ZZX b1;
      negate(b1, b);
      PseudoDivRem(q, r, a, b1);
      negate(q, q);
   }
   else if (divide(q, a, b)) {
      r = 0;
   }
   else {
      ZZX q1, r1;
      ZZ m;
      PseudoDivRem(q1, r1, a, b);
      power(m, LeadCoeff(b), da-db+1);
      if (!divide(q, q1, m)) Error("DivRem: quotient not defined over ZZ");
      if (!divide(r, r1, m)) Error("DivRem: remainder not defined over ZZ");
   }
}

void div(ZZX& q, const ZZX& a, const ZZX& b)
{
   long da = deg(a);
   long db = deg(b);

   if (db < 0) Error("div: division by zero");

   if (da < db) {
      q = 0;
   }
   else if (db == 0) {
      div(q, a, ConstTerm(b));
   }
   else if (IsOne(LeadCoeff(b))) {
      PseudoDiv(q, a, b);
   }
   else if (LeadCoeff(b) == -1) {
      ZZX b1;
      negate(b1, b);
      PseudoDiv(q, a, b1);
      negate(q, q);
   }
   else if (divide(q, a, b)) {

      // nothing to do
      
   }
   else {
      ZZX q1;
      ZZ m;
      PseudoDiv(q1, a, b);
      power(m, LeadCoeff(b), da-db+1);
      if (!divide(q, q1, m)) Error("div: quotient not defined over ZZ");
   }
}

void rem(ZZX& r, const ZZX& a, const ZZX& b)
{
   long da = deg(a);
   long db = deg(b);

   if (db < 0) Error("rem: division by zero");

   if (da < db) {
      r = a;
   }
   else if (db == 0) {
      ConstRem(r, a, ConstTerm(b));
   }
   else if (IsOne(LeadCoeff(b))) {
      PseudoRem(r, a, b);
   }
   else if (LeadCoeff(b) == -1) {
      ZZX b1;
      negate(b1, b);
      PseudoRem(r, a, b1);
   }
   else if (divide(a, b)) {
      r = 0;
   }
   else {
      ZZX r1;
      ZZ m;
      PseudoRem(r1, a, b);
      power(m, LeadCoeff(b), da-db+1);
      if (!divide(r, r1, m)) Error("rem: remainder not defined over ZZ");
   }
}

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

   if (IsZero(a)) {
      clear(q);
      return 1;
   }

   if (deg(b) == 0) {
      return divide(q, a, ConstTerm(b));
   }

   if (deg(a) < deg(b)) return 0;

   ZZ ca, cb, cq;

   content(ca, a);
   content(cb, b);

   if (!divide(cq, ca, cb)) return 0;

   ZZX aa, bb;

   divide(aa, a, ca);
   divide(bb, b, cb);

   if (!divide(LeadCoeff(aa), LeadCoeff(bb)))
      return 0;

   if (!divide(ConstTerm(aa), ConstTerm(bb)))
      return 0;

   zz_pBak bak;
   bak.save();

   ZZX qq;

   ZZ prod;
   set(prod);

   clear(qq);
   long res = 1;
   long Qinstable = 1;


   long a_bound = MaxBits(aa);
   long b_bound = MaxBits(bb);


   long i;
   for (i = 0; ; i++) {
      zz_p::FFTInit(i);
      long p = zz_p::modulus();

      if (divide(LeadCoeff(bb), p)) continue;

      zz_pX A, B, Q, R;

      conv(A, aa);
      conv(B, bb);

      if (!Qinstable) {
         conv(Q, qq);
         mul(R, B, Q);
         sub(R, A, R);

         if (deg(R) >= deg(B))
            Qinstable = 1;
         else if (!IsZero(R)) {
            res = 0;
            break;
         }
         else
            mul(prod, prod, p);
      }

      if (Qinstable) {
         if (!divide(Q, A, B)) {
            res = 0;
            break;
         }

         Qinstable = CRT(qq, prod, Q);
      }

      if (!Qinstable) {
         // stabilized...check if prod is big enough

         long bound = b_bound + MaxBits(qq) + 
                     NumBits(min(deg(bb), deg(qq)) + 1);

         if (a_bound > bound)
            bound = a_bound;

         bound += 3;

         if (NumBits(prod) > bound) 
            break;
      }
   }

   bak.restore();

   if (res) mul(q, qq, cq);
   return res;

}


long HomDivide(const ZZX& a, const ZZX& b)
{
   if (deg(b) == 0) {
      return divide(a, ConstTerm(b));
   }
   else {
      ZZX q;
      return HomDivide(q, a, b);
   }
}

long PlainDivide(ZZX& qq, const ZZX& aa, const ZZX& bb)
{
   if (IsZero(bb)) {
      if (IsZero(aa)) {
         clear(qq);
         return 1;
      }
      else
         return 0;
   }

   if (deg(bb) == 0) {
      return divide(qq, aa, ConstTerm(bb));
   }

   long da, db, dq, i, j, LCIsOne;
   const ZZ *bp;
   ZZ *qp;
   ZZ *xp;


   ZZ  s, t;

   da = deg(aa);
   db = deg(bb);

   if (da < db) {
      return 0;
   }

   ZZ ca, cb, cq;

   content(ca, aa);
   content(cb, bb);

   if (!divide(cq, ca, cb)) {
      return 0;
   } 


   ZZX a, b, q;

   divide(a, aa, ca);
   divide(b, bb, cb);

   if (!divide(LeadCoeff(a), LeadCoeff(b)))
      return 0;

   if (!divide(ConstTerm(a), ConstTerm(b)))
      return 0;

   long coeff_bnd = MaxBits(a) + (NumBits(da+1)+1)/2 + (da-db);

   bp = b.rep.elts();

   ZZ LC;
   LC = bp[db];

   LCIsOne = IsOne(LC);

   xp = a.rep.elts();

   dq = da - db;
   q.rep.SetLength(dq+1);
   qp = q.rep.elts();

   for (i = dq; i >= 0; i--) {
      if (!LCIsOne) {
         if (!divide(t, xp[i+db], LC))
            return 0;
      }
      else
         t = xp[i+db];

      if (NumBits(t) > coeff_bnd) return 0;

      qp[i] = t;

      for (j = db-1; j >= 0; j--) {
	 mul(s, t, bp[j]);
	 sub(xp[i+j], xp[i+j], s);
      }
   }

   for (i = 0; i < db; i++)
      if (!IsZero(xp[i]))
         return 0;

   mul(qq, q, cq);
   return 1;
}

long PlainDivide(const ZZX& a, const ZZX& b)
{
   if (deg(b) == 0) 
      return divide(a, ConstTerm(b));
   else {
      ZZX q;
      return PlainDivide(q, a, b);
   }
}


long divide(ZZX& q, const ZZX& a, const ZZX& b)
{
   long da = deg(a);
   long db = deg(b);

   if (db <= 8 || da-db <= 8)
      return PlainDivide(q, a, b);
   else
      return HomDivide(q, a, b);
}

long divide(const ZZX& a, const ZZX& b)
{
   long da = deg(a);
   long db = deg(b);

   if (db <= 8 || da-db <= 8)
      return PlainDivide(a, b);
   else
      return HomDivide(a, b);
}







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

   if (IsOne(b)) {
      q = a;
      return 1;
   }

   if (b == -1) {
      negate(q, a);
      return 1;
   }

   long n = a.rep.length();
   vec_ZZ res(INIT_SIZE, n);
   long i;

   for (i = 0; i < n; i++) {
      if (!divide(res[i], a.rep[i], b))
         return 0;
   }

   q.rep = res;
   return 1;
}

long divide(const ZZX& a, const ZZ& b)
{
   if (IsZero(b)) return IsZero(a);

   if (IsOne(b) || b == -1) {
      return 1;
   }

   long n = a.rep.length();
   long i;

   for (i = 0; i < n; i++) {
      if (!divide(a.rep[i], b))
         return 0;
   }

   return 1;
}

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

   if (b == 1) {
      q = a;
      return 1;
   }

   if (b == -1) {
      negate(q, a);
      return 1;
   }

   long n = a.rep.length();
   vec_ZZ res(INIT_SIZE, n);
   long i;

   for (i = 0; i < n; i++) {
      if (!divide(res[i], a.rep[i], b))
         return 0;
   }

   q.rep = res;
   return 1;
}

long divide(const ZZX& a, long b)
{
   if (b == 0) return IsZero(a);
   if (b == 1 || b == -1) {
      return 1;
   }

   long n = a.rep.length();
   long i;

   for (i = 0; i < n; i++) {
      if (!divide(a.rep[i], b))
         return 0;
   }

   return 1;
}

   

void content(ZZ& d, const ZZX& f)
{
   ZZ res;
   long i;

   clear(res);
   for (i = 0; i <= deg(f); i++) {
      GCD(res, res, f.rep[i]);
      if (IsOne(res)) break;
   }

   if (sign(LeadCoeff(f)) < 0) negate(res, res);
   d = res;
}

void PrimitivePart(ZZX& pp, const ZZX& f)
{
   if (IsZero(f)) {
      clear(pp);
      return;
   }
 
   ZZ d;

   content(d, f);
   divide(pp, f, d);
}


static
void BalCopy(ZZX& g, const zz_pX& G)
{
   long p = zz_p::modulus();
   long p2 = p >> 1;
   long n = G.rep.length();
   long i;
   long t;

   g.rep.SetLength(n);
   for (i = 0; i < n; i++) {
      t = rep(G.rep[i]);
      if (t > p2) t = t - p;
      conv(g.rep[i], t);
   }
}


   
void GCD(ZZX& d, const ZZX& a, const ZZX& b)
{
   if (IsZero(a)) {
      d = b;
      if (sign(LeadCoeff(d)) < 0) negate(d, d);
      return;
   }

   if (IsZero(b)) {
      d = a;
      if (sign(LeadCoeff(d)) < 0) negate(d, d);
      return;
   }

   ZZ c1, c2, c;
   ZZX f1, f2;

   content(c1, a);
   divide(f1, a, c1);

   content(c2, b);
   divide(f2, b, c2);

   GCD(c, c1, c2);

   ZZ ld;
   GCD(ld, LeadCoeff(f1), LeadCoeff(f2));

   ZZX g, h, res;