poly.tcc

CGAL is a collaborative effort of several sites in Europe and Israel. The goal is to make the most i

  if(p.getTrueDegree() < q.getTrueDegree())
    return false;
  p.pseudoRemainder(q);
  if(zeroP(p))
    return true;
  else
    return false;
}//isDivisible

//Content of a polynomial P
//      -- content(P) is just the gcd of all the coefficients
//      -- REMARK: by definition, content(P) is non-negative
//                 We rely on the fact that NT::gcd always
//                 return a non-negative value!!!
template <class NT>
NT content(const Polynomial<NT>& p) {
  if(zeroP(p))
    return 0;
  int d = p.getTrueDegree();
  if(d == 0){
    if(p.coeff[0] > 0)
      return p.coeff[0];
    else
      return -p.coeff[0];
  }

  NT content = p.coeff[d];
  for (int i=d-1; i>=0; i--) {
    content = gcd(content, p.coeff[i]);
    if(content == 1) break;   // remark: gcd is non-negative, by definition
  }
  //if (p.coeff[0] < 0) return -content;(BUG!)
  return content;
}//content

// Primitive Part:  (*this) is transformed to primPart and returned
//	-- primPart(P) is just P/content(P)
//	-- Should we return content(P) instead? [SHOULD IMPLEMENT THIS]
// IMPORTANT: we require that content(P)>0, hence
// 	the coefficients of primPart(P) does 
// 	not change sign; this is vital for use in Sturm sequences
template <class NT>
Polynomial<NT> & Polynomial<NT>::primPart() {
  // ASSERT: GCD must be provided by NT
  int d = getTrueDegree();
  assert (d >= 0);
  if (d == 0) {
    if (coeff[0] > 0) coeff[0] = 1;
    else coeff[0] = -1;
    return *this;
  }

  NT g = content(*this);
  if (g == 1 && coeff[d] > 0)
     return (*this);  
  for (int i=0; i<=d; i++) {
     coeff[i] =  coeff[i]/g;
  }
  return *this;
}// primPart

//GCD of two polynomials.
//   --Assumes that the coeffecient ring NT has a gcd function
//   --Returns the gcd with a positive leading coefficient (*)
//     otherwise division by gcd causes a change of sign affecting
//     Sturm sequences.
//   --To Check: would a non-recursive version be much better?
template <class NT>
Polynomial<NT> gcd(const Polynomial<NT>& p, const Polynomial<NT>& q) {

  // If the first polynomial has a smaller degree then the second,
  // then change the order of calling
  if(p.getTrueDegree() < q.getTrueDegree())
    return gcd(q,p);

  // If any polynomial is zero then the gcd is the other polynomial
  if(zeroP(q)){
    if(zeroP(p))
       return p;
    else{
       if(p.getCoeffi(p.getTrueDegree()) < 0){
         return Polynomial<NT>(p).negate();
       }else
         return p;	// If q<>0, then we know p<>0
   }
  }
  Polynomial<NT> temp0(p);
  Polynomial<NT> temp1(q);

  // We want to compute:
  // gcd(p,q) = gcd(content(p),content(q)) * gcd(primPart(p), primPart(q))

  NT cont0 = content(p);	// why is this temporary needed?
  NT cont1 = content(q);
  NT cont = gcd(cont0,cont1);
  temp0.primPart();
  temp1.primPart();

  temp0.pseudoRemainder(temp1);
  return (gcd(temp1, temp0).mulScalar(cont));
}//gcd

// sqFreePart()
// 	-- this is a self-modifying operator!
// 	-- Let P =(*this) and Q=square-free part of P.
// 	-- (*this) is transformed into P, and gcd(P,P') is returned
// NOTE: The square-free part of P is defined to be P/gcd(P, P')
template <class NT>
Polynomial<NT>  Polynomial<NT>::sqFreePart() {

  int d = getTrueDegree();
  if(d <= 1) // linear polynomials and constants are square-free
    return *this;

  Polynomial<NT> temp(*this);
  Polynomial<NT> R = gcd(*this, temp.differentiate()); // R = gcd(P, P')

  // If P and P' have a constant gcd, then P is squarefree
  if(R.getTrueDegree() == 0)
    return (Polynomial<NT>(0)); // returns the unit polynomial as gcd

  (*this)=pseudoRemainder(R); // (*this) is transformed to P/R, the sqfree part
  //Note: This is up to multiplication by units
  return (R); // return the gcd
}//sqFreePart()


// ==================================================
// Useful member functions
// ==================================================

// reverse:
// 	reverses the list of coefficients
template <class NT>
void Polynomial<NT>::reverse() {
  NT tmp;
  for (int i=0; i<= degree/2; i++) {
    tmp = coeff[i];
    coeff[i] =   coeff[degree-i];
    coeff[degree-i] = tmp;
  }
}//reverse

// negate: 
// 	multiplies the polynomial by -1
// 	Chee: 4/29/04 -- added negate() to support negPseudoRemainder(B)
template <class NT>
Polynomial<NT> & Polynomial<NT>::negate() {
  for (int i=0; i<= degree; i++) 
    coeff[i] *= -1;  	// every NT must be able to construct from -1
  return *this;
}//negate

// makeTailCoeffNonzero
// 	Divide (*this) by X^k, so that the tail coeff becomes non-zero.
//	The value k is returned.  In case (*this) is 0, return k=-1.
//	Otherwise, if (*this) is unchanged, return k=0.
template <class NT>
int Polynomial<NT>::makeTailCoeffNonzero() {
  int k=-1;
  for (int i=0; i <= degree; i++) {
    if (coeff[i] != 0) {
      k=i;
      break;
    }
  }
  if (k <= 0)
    return k;	// return either 0 or -1
  degree -=k;		// new (lowered) degree
  NT * c = new NT[1+degree];
  for (int i=0; i<=degree; i++)
    c[i] = coeff[i+k];
  delete[] coeff;
  coeff = c;
  return k;
}//

// filedump(string msg, ostream os, string com, string com2):
//      Dumps polynomial to output stream os
//      msg is any message
//      NOTE: Default is com="#", which is placed at start of each 
//            output line. 
template <class NT>
void Polynomial<NT>::filedump(std::ostream & os,
                          std::string msg,
			  std::string commentString,
                          std::string commentString2) const {
  int d= getTrueDegree();
  if (msg != "") os << commentString << msg << std::endl;
  int i=0;
  if (d == -1) { // if zero polynomial
    os << commentString << "0";
    return;
  }
  for (; i<=d;  ++i)	// get to the first non-zero coeff
    if (coeff[i] != 0)
      break;
  int termsInLine = 1;

  // OUTPUT the first nonzero term
  os << commentString;
  if (coeff[i] == 1) {			// special cases when coeff[i] is
    if (i>1) os << "x^" << i;	// either 1 or -1 
    else if (i==1) os << "x" ;
    else os << "1";
  } else if (coeff[i] == -1) {
    if (i>1) os << "-x^" << i;
    else if (i==1) os << "-x" ;
    else os << "-1";
  } else {				// non-zero coeff
    os << coeff[i];
    if (i>1) os << "*x^" << i;
    else if (i==1) os << "x" ;
  } 
  // OUTPUT the remaining nonzero terms
  for (i++ ; i<= getTrueDegree(); ++i) {
    if (coeff[i] == 0) 
      continue;
    termsInLine++;
    if (termsInLine % Polynomial<NT>::COEFF_PER_LINE == 0) {
      os << std::endl;
      os << commentString2;
    }
    if (coeff[i] == 1) {		// special when coeff[i] = 1
      if (i==1) os << " + x";
      else os << " + x^" << i;
    } else if (coeff[i] == -1) {	// special when coeff[i] = -1
      if (i==1) os << " - x";
      else os << " -x^" << i;
    } else {				// general case
      if(coeff[i] > 0){	
	os << " + ";
        os << coeff[i];
      }else
        os << coeff[i];

      if (i==1) os << "*x";
      else os << "*x^" << i;
    } 
  }
}//filedump

// dump(message, ofstream, commentString) -- dump to file
template <class NT>
void Polynomial<NT>::dump(std::ofstream & ofs,
		std::string msg,
		std::string commentString,
                std::string commentString2) const {
  filedump(ofs, msg, commentString, commentString2);
}

// dump(message) 	-- to std output
template <class NT>
void Polynomial<NT>::dump(std::string msg, std::string com,
		std::string com2) const {
  filedump(std::cout, msg, com, com2);
}

// Dump of Maple Code for Polynomial
template <class NT>
void Polynomial<NT>::mapleDump() const {
  if (zeroP(*this)) {
    std::cout << 0 << std::endl;
    return;
  }
  std::cout << coeff[0];
  for (int i = 1; i<= getTrueDegree(); ++i) {
    std::cout << " + (" << coeff[i] << ")";
    std::cout << "*x^" << i;
  }
  std::cout << std::endl;
}//mapleDump

// ==================================================
// Useful friend functions for Polynomial<NT> class
// ==================================================

// friend differentiation
template <class NT>
Polynomial<NT> differentiate(const Polynomial<NT> & p) {	  // differentiate
  Polynomial<NT> q(p);
  return q.differentiate();
}

// friend multi-differentiation
template <class NT>
Polynomial<NT> differentiate(const Polynomial<NT> & p, int n) {//multi-differentiate
  Polynomial<NT> q(p);
  assert(n >= 0);
  for (int i=1; i<=n; i++)
    q.differentiate();
  return q;
}

// friend equality comparison
template <class NT>
bool operator==(const Polynomial<NT>& p, const Polynomial<NT>& q) {	// ==
  int d, D;
  Polynomial<NT> P(p);
  P.contract();
  Polynomial<NT> Q(q);
  Q.contract();
  if (P.degree < Q.degree) {
    d = P.degree;
    D = Q.degree;
    for (int i = d+1; i<=D; i++)
      if (Q.coeff[i] != 0)
        return false;	// return false
  } else {
    D = P.degree;
    d = Q.degree;
    for (int i = d+1; i<=D; i++)
      if (P.coeff[i] != 0)
        return false;	// return false
  }
  for (int i = 0; i <= d; i++)
    if (P.coeff[i] != Q.coeff[i])
      return false;	// return false
  return true; 	// return true
}

// friend non-equality comparison
template <class NT>
bool operator!=(const Polynomial<NT>& p, const Polynomial<NT>& q) {	// !=
  return (!(p == q));
}

// stream i/o
template <class NT>
std::ostream& operator<<(std::ostream& o, const Polynomial<NT>& p) {
  o <<   "Polynomial<NT> ( deg = " << p.degree ;
  if (p.degree >= 0) {
    o << "," << std::endl;
    o << ">  coeff c0,c1,... = " << p.coeff[0];
    for (int i=1; i<= p.degree; i++)
      o << ", " <<  p.coeff[i] ;
  }
  o << ")" << std::endl;
  return o;
}

// fragile version...
template <class NT>
std::istream& operator>>(std::istream& is, Polynomial<NT>& p) {
  is >> p.degree;
  // Don't you need to first do  "delete[] p.coeff;"  ??
  p.coeff = new NT[p.degree+1];
  for (int i=0; i<= p.degree; i++)
    is >> p.coeff[i];
  return is;
}

// ==================================================
// Simple test of poly
// ==================================================

template <class NT>
bool testPoly() {
  int c[] = {1, 2, 3};
  Polynomial<NT> p(2, c);
  std::cout << p;

  Polynomial<NT> zeroP;
  std::cout << "zeroP  : " << zeroP << std::endl;

  Polynomial<NT> P5(5);
  std::cout << "Poly 5 : " << P5 << std::endl;

  return 0;
}


//Resultant of two polynomials.
//Since we use pseudoRemainder we have to modify the original algorithm.
//If C * P = Q * R + S, where C is a constant and S = prem(P,Q), m=deg(P),
// n=deg(Q) and l = deg(S).
//Then res(P,Q) = (-1)^(mn) b^(m-l) res(Q,S)/C^(n)
//The base case being res(P,Q) = Q^(deg(P)) if Q is a constant, zero otherwise
template <class NT>
NT res(Polynomial<NT> p, Polynomial<NT> q) {

  int m, n;
  m = p.getTrueDegree();
  n = q.getTrueDegree();

  if(m == -1 || n == -1) return 0;  // this definition is not certified
  if(m == 0 && n == 0) return 1;    // this definition is at variance from
                                    // Yap's book (but is OK for purposes of
                                    // checking the vanishing of resultants
  if(n > m) return (res(q, p));

  NT b = q.coeff[n];//Leading coefficient of Q
  NT lc = p.coeff[m], C;

  p.pseudoRemainder(q, C);

  if(zeroP(p) && n ==0)
     return (pow(q.coeff[0], m));	

  int l = p.getTrueDegree();

  return(pow(NT(-1), m*n)*pow(b,m-l)*res(q,p)/pow(C,n));
}

//i^th Principal Subresultant Coefficient (psc) of two polynomials.
template <class NT>
NT psc(int i,Polynomial<NT> p, Polynomial<NT> q) {

  assert(i >= 0);
  if(i == 0)
     return res(p,q);

  int m = p.getTrueDegree();
  int n = q.getTrueDegree();

  if(m == -1 || n == -1) return 0;
  if(m < n) return psc(i, q, p);

  if ( i == n) //This case occurs when i > degree of gcd
     return pow(q.coeff[n], m - n);

  if(n < i && i <= m)
     return 0;

  NT b = q.coeff[n];//Leading coefficient of Q
  NT lc = p.coeff[m], C;	


  p.pseudoRemainder(q, C);

  if(zeroP(p) && i < n)//This case occurs when i < deg(gcd)
     return 0;

  if(zeroP(p) && i == n)//This case occurs when i=deg(gcd) might be constant
     return pow(q.coeff[n], m - n);

  int l = p.getTrueDegree();
  return pow(NT(-1),(m-i)*(n-i))*pow(b,m-l)*psc(i,q,p)/pow(C, n-i);

}

//factorI(p,m)
//    computes the polynomial q which containing all roots
//    of multiplicity m of a given polynomial P
//    Used to determine the nature of intersection of two algebraic curves
//    The algorithm is given in Wolperts Thesis
template <class NT>
Polynomial<NT> factorI(Polynomial<NT> p, int m){
  int d=p.getTrueDegree();
  Polynomial<NT> *w = new Polynomial<NT>[d+1];
  w[0] = p;
  Polynomial<NT> temp;

  for(int i = 1; i <=d ; i ++){
     temp = w[i-1];
     w[i] = gcd(w[i-1],temp.differentiate());
  }

  Polynomial<NT> *u = new Polynomial<NT>[d+1];
  u[d] = w[d-1];
  for(int i = d-1; i >=m; i--){
	temp = power(u[i+1],2);
     for(int j=i+2; j <=d; j++){
        temp *= power(u[j],j-i+1);
     }
     u[i] = w[i-1].pseudoRemainder(temp);//i-1 since the array starts at 0
  } 
	
	delete[] w;
  return u[m];
}

// ==================================================
// End of Polynomial<NT>
// ==================================================