lzz_pex.cpp

数值算法库for Windows

#include <NTL/lzz_pEX.h>
#include <NTL/vec_vec_lzz_p.h>

#include <NTL/new.h>

NTL_START_IMPL


const zz_pEX& zz_pEX::zero()
{
   static zz_pEX z;
   return z;
}


istream& operator>>(istream& s, zz_pEX& x)
{
   s >> x.rep;
   x.normalize();
   return s;
}

ostream& operator<<(ostream& s, const zz_pEX& a)
{
   return s << a.rep;
}


void zz_pEX::normalize()
{
   long n;
   const zz_pE* p;

   n = rep.length();
   if (n == 0) return;
   p = rep.elts() + (n-1);
   while (n > 0 && IsZero(*p)) {
      p--; 
      n--;
   }
   rep.SetLength(n);
}


long IsZero(const zz_pEX& a)
{
   return a.rep.length() == 0;
}


long IsOne(const zz_pEX& a)
{
    return a.rep.length() == 1 && IsOne(a.rep[0]);
}

long operator==(const zz_pEX& a, long b)
{
   if (b == 0)
      return IsZero(a);

   if (b == 1)
      return IsOne(a);

   long da = deg(a);

   if (da > 0) return 0;

   NTL_zz_pRegister(bb);
   bb = b;

   if (da < 0)
      return IsZero(bb);

   return a.rep[0] == bb;
}

long operator==(const zz_pEX& a, const zz_p& b)
{
   if (IsZero(b))
      return IsZero(a);

   long da = deg(a);

   if (da != 0)
      return 0;

   return a.rep[0] == b;
}

long operator==(const zz_pEX& a, const zz_pE& b)
{
   if (IsZero(b))
      return IsZero(a);

   long da = deg(a);

   if (da != 0)
      return 0;

   return a.rep[0] == b;
}




void SetCoeff(zz_pEX& x, long i, const zz_pE& a)
{
   long j, m;

   if (i < 0) 
      Error("SetCoeff: negative index");

   if (i >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("overflow in SetCoeff");

   m = deg(x);

   if (i > m) {
      long pos = x.rep.position(a);
      x.rep.SetLength(i+1);

      if (pos != -1)
         x.rep[i] = x.rep.RawGet(pos);
      else
         x.rep[i] = a;

      for (j = m+1; j < i; j++)
         clear(x.rep[j]);
   }
   else
      x.rep[i] = a;

   x.normalize();
}

void SetCoeff(zz_pEX& x, long i, const zz_p& aa)
{
   long j, m;

   if (i < 0)
      Error("SetCoeff: negative index");

   if (i >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("overflow in SetCoeff");

   NTL_zz_pRegister(a);  // watch out for aliases!
   a = aa;

   m = deg(x);

   if (i > m) {
      x.rep.SetLength(i+1);
      for (j = m+1; j < i; j++)
         clear(x.rep[j]);
   }
   x.rep[i] = a;
   x.normalize();
}

void SetCoeff(zz_pEX& x, long i, long a)
{
   if (a == 1)
      SetCoeff(x, i);
   else {
      NTL_zz_pRegister(T);
      T = a;
      SetCoeff(x, i, T);
   }
}



void SetCoeff(zz_pEX& x, long i)
{
   long j, m;

   if (i < 0) 
      Error("coefficient index out of range");

   if (i >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("overflow in SetCoeff");

   m = deg(x);

   if (i > m) {
      x.rep.SetLength(i+1);
      for (j = m+1; j < i; j++)
         clear(x.rep[j]);
   }
   set(x.rep[i]);
   x.normalize();
}


void SetX(zz_pEX& x)
{
   clear(x);
   SetCoeff(x, 1);
}


long IsX(const zz_pEX& a)
{
   return deg(a) == 1 && IsOne(LeadCoeff(a)) && IsZero(ConstTerm(a));
}
      
      

const zz_pE& coeff(const zz_pEX& a, long i)
{
   if (i < 0 || i > deg(a))
      return zz_pE::zero();
   else
      return a.rep[i];
}


const zz_pE& LeadCoeff(const zz_pEX& a)
{
   if (IsZero(a))
      return zz_pE::zero();
   else
      return a.rep[deg(a)];
}

const zz_pE& ConstTerm(const zz_pEX& a)
{
   if (IsZero(a))
      return zz_pE::zero();
   else
      return a.rep[0];
}



void conv(zz_pEX& x, const zz_pE& a)
{
   if (IsZero(a))
      x.rep.SetLength(0);
   else {
      x.rep.SetLength(1);
      x.rep[0] = a;
   }
}

void conv(zz_pEX& x, long a)
{
   if (a == 0) 
      clear(x);
   else if (a == 1)
      set(x);
   else {
      NTL_zz_pRegister(T);
      T = a;
      conv(x, T);
   }
}

void conv(zz_pEX& x, const ZZ& a)
{
   NTL_zz_pRegister(T);
   conv(T, a);
   conv(x, T);
}

void conv(zz_pEX& x, const zz_p& a)
{
   if (IsZero(a)) 
      clear(x);
   else if (IsOne(a))
      set(x);
   else {
      x.rep.SetLength(1);
      conv(x.rep[0], a);
      x.normalize();
   }
}

void conv(zz_pEX& x, const zz_pX& aa)
{
   zz_pX a = aa; // in case a aliases the rep of a coefficient of x

   long n = deg(a)+1;
   long i;

   x.rep.SetLength(n);
   for (i = 0; i < n; i++)
      conv(x.rep[i], coeff(a, i));
}


void conv(zz_pEX& x, const vec_zz_pE& a)
{
   x.rep = a;
   x.normalize();
}


void add(zz_pEX& x, const zz_pEX& a, const zz_pEX& b)
{
   long da = deg(a);
   long db = deg(b);
   long minab = min(da, db);
   long maxab = max(da, db);
   x.rep.SetLength(maxab+1);

   long i;
   const zz_pE *ap, *bp; 
   zz_pE* xp;

   for (i = minab+1, ap = a.rep.elts(), bp = b.rep.elts(), xp = x.rep.elts();
        i; i--, ap++, bp++, xp++)
      add(*xp, (*ap), (*bp));

   if (da > minab && &x != &a)
      for (i = da-minab; i; i--, xp++, ap++)
         *xp = *ap;
   else if (db > minab && &x != &b)
      for (i = db-minab; i; i--, xp++, bp++)
         *xp = *bp;
   else
      x.normalize();
}


void add(zz_pEX& x, const zz_pEX& a, const zz_pE& b)
{
   long n = a.rep.length();
   if (n == 0) {
      conv(x, b);
   }
   else if (&x == &a) {
      add(x.rep[0], a.rep[0], b);
      x.normalize();
   }
   else if (x.rep.MaxLength() == 0) {
      x = a;
      add(x.rep[0], a.rep[0], b);
      x.normalize();
   }
   else {
      // ugly...b could alias a coeff of x

      zz_pE *xp = x.rep.elts();
      add(xp[0], a.rep[0], b);
      x.rep.SetLength(n);
      xp = x.rep.elts();
      const zz_pE *ap = a.rep.elts();
      long i;
      for (i = 1; i < n; i++)
         xp[i] = ap[i];
      x.normalize();
   }
}

void add(zz_pEX& x, const zz_pEX& a, const zz_p& b)
{
   long n = a.rep.length();
   if (n == 0) {
      conv(x, b);
   }
   else if (&x == &a) {
      add(x.rep[0], a.rep[0], b);
      x.normalize();
   }
   else if (x.rep.MaxLength() == 0) {
      x = a;
      add(x.rep[0], a.rep[0], b);
      x.normalize();
   }
   else {
      // ugly...b could alias a coeff of x

      zz_pE *xp = x.rep.elts();
      add(xp[0], a.rep[0], b);
      x.rep.SetLength(n);
      xp = x.rep.elts();
      const zz_pE *ap = a.rep.elts();
      long i;
      for (i = 1; i < n; i++)
         xp[i] = ap[i];
      x.normalize();
   }
}


void add(zz_pEX& x, const zz_pEX& a, long b)
{
   if (a.rep.length() == 0) {
      conv(x, b);
   }
   else {
      if (&x != &a) x = a;
      add(x.rep[0], x.rep[0], b);
      x.normalize();
   }
}


void sub(zz_pEX& x, const zz_pEX& a, const zz_pEX& b)
{
   long da = deg(a);
   long db = deg(b);
   long minab = min(da, db);
   long maxab = max(da, db);
   x.rep.SetLength(maxab+1);

   long i;
   const zz_pE *ap, *bp; 
   zz_pE* xp;

   for (i = minab+1, ap = a.rep.elts(), bp = b.rep.elts(), xp = x.rep.elts();
        i; i--, ap++, bp++, xp++)
      sub(*xp, (*ap), (*bp));

   if (da > minab && &x != &a)
      for (i = da-minab; i; i--, xp++, ap++)
         *xp = *ap;
   else if (db > minab)
      for (i = db-minab; i; i--, xp++, bp++)
         negate(*xp, *bp);
   else
      x.normalize();
}


void sub(zz_pEX& x, const zz_pEX& a, const zz_pE& b)
{
   long n = a.rep.length();
   if (n == 0) {
      conv(x, b);
      negate(x, x);
   }
   else if (&x == &a) {
      sub(x.rep[0], a.rep[0], b);
      x.normalize();
   }
   else if (x.rep.MaxLength() == 0) {
      x = a;
      sub(x.rep[0], a.rep[0], b);
      x.normalize();
   }
   else {
      // ugly...b could alias a coeff of x

      zz_pE *xp = x.rep.elts();
      sub(xp[0], a.rep[0], b);
      x.rep.SetLength(n);
      xp = x.rep.elts();
      const zz_pE *ap = a.rep.elts();
      long i;
      for (i = 1; i < n; i++)
         xp[i] = ap[i];
      x.normalize();
   }
}

void sub(zz_pEX& x, const zz_pEX& a, const zz_p& b)
{
   long n = a.rep.length();
   if (n == 0) {
      conv(x, b);
      negate(x, x);
   }
   else if (&x == &a) {
      sub(x.rep[0], a.rep[0], b);
      x.normalize();
   }
   else if (x.rep.MaxLength() == 0) {
      x = a;
      sub(x.rep[0], a.rep[0], b);
      x.normalize();
   }
   else {
      // ugly...b could alias a coeff of x

      zz_pE *xp = x.rep.elts();
      sub(xp[0], a.rep[0], b);
      x.rep.SetLength(n);
      xp = x.rep.elts();
      const zz_pE *ap = a.rep.elts();
      long i;
      for (i = 1; i < n; i++)
         xp[i] = ap[i];
      x.normalize();
   }
}


void sub(zz_pEX& x, const zz_pEX& a, long b)
{
   if (a.rep.length() == 0) {
      conv(x, b);
      negate(x, x);
   }
   else {
      if (&x != &a) x = a;
      sub(x.rep[0], x.rep[0], b);
      x.normalize();
   }
}

void sub(zz_pEX& x, const zz_pE& b, const zz_pEX& a)
{
   long n = a.rep.length();
   if (n == 0) {
      conv(x, b);
   }
   else if (x.rep.MaxLength() == 0) {
      negate(x, a);
      add(x.rep[0], a.rep[0], b);
      x.normalize();
   }
   else {
      // ugly...b could alias a coeff of x

      zz_pE *xp = x.rep.elts();
      sub(xp[0], b, a.rep[0]);
      x.rep.SetLength(n);
      xp = x.rep.elts();
      const zz_pE *ap = a.rep.elts();
      long i;
      for (i = 1; i < n; i++)
         negate(xp[i], ap[i]);
      x.normalize();
   }
}


void sub(zz_pEX& x, const zz_p& a, const zz_pEX& b)
{
   NTL_zz_pRegister(T);   // avoids aliasing problems
   T = a;
   negate(x, b);
   add(x, x, T);
}

void sub(zz_pEX& x, long a, const zz_pEX& b)
{
   NTL_zz_pRegister(T); 
   T = a;
   negate(x, b);
   add(x, x, T);
}

void mul(zz_pEX& c, const zz_pEX& a, const zz_pEX& b)
{
   if (&a == &b) {
      sqr(c, a);
      return;
   }

   if (IsZero(a) || IsZero(b)) {
      clear(c);
      return;
   }

   if (deg(a) == 0) {
      mul(c, b, ConstTerm(a));
      return;
   } 

   if (deg(b) == 0) {
      mul(c, a, ConstTerm(b));
      return;
   }

   // general case...Kronecker subst

   zz_pX A, B, C;

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

   long n = zz_pE::degree();
   long n2 = 2*n-1;

   if (da+db+1 >= (1L << (NTL_BITS_PER_LONG-4))/n2)
      Error("overflow in zz_pEX mul");


   long i, j;

   A.rep.SetLength((da+1)*n2);

   for (i = 0; i <= da; i++) {
      const zz_pX& coeff = rep(a.rep[i]);
      long dcoeff = deg(coeff);
      for (j = 0; j <= dcoeff; j++)
         A.rep[n2*i + j] = coeff.rep[j]; 
   }

   A.normalize();

   B.rep.SetLength((db+1)*n2);

   for (i = 0; i <= db; i++) {
      const zz_pX& coeff = rep(b.rep[i]);
      long dcoeff = deg(coeff);
      for (j = 0; j <= dcoeff; j++)
         B.rep[n2*i + j] = coeff.rep[j]; 
   }

   B.normalize();

   mul(C, A, B);

   long Clen = C.rep.length();
   long lc = (Clen + n2 - 1)/n2;
   long dc = lc - 1;

   c.rep.SetLength(dc+1);

   zz_pX tmp;
   
   for (i = 0; i <= dc; i++) {
      tmp.rep.SetLength(n2);
      for (j = 0; j < n2 && n2*i + j < Clen; j++)
         tmp.rep[j] = C.rep[n2*i + j];
      for (; j < n2; j++)
         clear(tmp.rep[j]);
      tmp.normalize();
      conv(c.rep[i], tmp);
   }
  
   c.normalize();
}


void mul(zz_pEX& x, const zz_pEX& a, const zz_pE& b)
{
   if (IsZero(b)) {
      clear(x);
      return;
   }

   zz_pE t;
   t = b;

   long i, da;

   const zz_pE *ap;
   zz_pE* xp;

   da = deg(a);
   x.rep.SetLength(da+1);
   ap = a.rep.elts();
   xp = x.rep.elts();

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

   x.normalize();
}



void mul(zz_pEX& x, const zz_pEX& a, const zz_p& b)
{
   if (IsZero(b)) {
      clear(x);
      return;
   }

   NTL_zz_pRegister(t);
   t = b;

   long i, da;

   const zz_pE *ap;
   zz_pE* xp;

   da = deg(a);
   x.rep.SetLength(da+1);
   ap = a.rep.elts();
   xp = x.rep.elts();