poly.tcc

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

// Multiply by a scalar
template <class NT>
Polynomial<NT> & Polynomial<NT>::mulScalar( const NT & c) {
  for (int i = 0; i<=degree ; i++)
    coeff[i] *= c;
  return *this;
}

// mulXpower: Multiply by X^i (COULD be a divide if i<0)
template <class NT>
Polynomial<NT> & Polynomial<NT>::mulXpower(int s) {
  // if s >= 0, then this is equivalent to
  // multiplying by X^s;  if s < 0, to dividing by X^s
  if (s==0)
    return *this;
  int d = s+getTrueDegree();
  if (d < 0) {
    degree = -1;
    delete[] coeff;
    coeff = NULL;
    return *this;
  }
  NT * c = new NT[d+1];
  if (s>0)
    for (int j=0;  j <= d; j++) {
      if (j <= degree)
        c[d-j] = coeff[d-s-j];
      else
        c[d-j] = 0;
    }
  if (s<0) {
    for (int j=0; j <= d; j++)
      c[d-j] = coeff[d-s-j];  // since s<0, (d-s-j) > (d-j) !!
  }
  delete[] coeff;
  coeff = c;
  degree = d;
  return *this;
}//mulXpower

// REDUCE STEP (helper for PSEUDO-REMAINDER function)
// Let THIS=(*this) be the current polynomial, and P be the input
//	argument for reduceStep.  Let R be returned polynomial.
//	R has the special form as a binomial,
//		R = C + X*M
//	where C is a constant and M a monomial (= coeff * some power of X).
//	Moreover, THIS is transformed to a new polynomial, THAT, which
//	is given by
// 		(C * THIS) = M * P  + THAT
//	MOREOVER: deg(THAT) < deg(THIS) unless deg(P)>deg(Q).
//	Basically, C is a power of the leading coefficient of P.
//	REMARK: R is NOT defined as C+M, because in case M is
//	a constant, then we cannot separate C from M.
//	Furthermore, R.mulScalar(-1) gives us M.
template <class NT>
Polynomial<NT> Polynomial<NT>::reduceStep (
  const Polynomial<NT>& p) {
  // 	Chee: Omit the next 2 contractions as unnecessary
  // 	since reduceStep() is only called by pseudoRemainder().
  // 	Also, reduceStep() already does a contraction before returning.
  // p.contract();	
  // contract();	// first contract both polynomials
  Polynomial<NT> q(p);		// q is initially a copy of p
  //	but is eventually M*P
  int pDeg  = q.degree;
  int myDeg = degree;
  if (pDeg == -1)
    return *(new Polynomial());  // Zero Polynomial
  // NOTE: pDeg=-1 (i.e., p=0) is really an error condition!
  if (myDeg < pDeg)
    return *(new Polynomial(0));  // Unity Polynomial
  // i.e., C=1, M=0.
  // Now (myDeg >= pDeg).  Start to form the Return Polynomial R=C+X*M
  Polynomial<NT> R(myDeg - pDeg + 1);  // deg(M)= myDeg - pDeg
  q.mulXpower(myDeg - pDeg);  	 // q is going to become M*P

  NT myLC = coeff[myDeg];	  // LC means "leading coefficient"
  NT qLC = q.coeff[myDeg];  // p also has degree "myDeg" (qLC non-zero)
  NT LC;

  //  NT must support
  //  isDivisible(x,y), gcd(x,y), div_exact(x,y) in the following:
  //  ============================================================
  if (isDivisible(myLC, qLC)) { // myLC is divisible by qLC
    LC = div_exact(myLC, qLC);	 
    R.setCoeff(0, 1);  		 //  C = 1,

    R.setCoeff(R.degree, LC); //  M = LC * X^(myDeg-pDeg)
    q.mulScalar(LC); 	  //  q = M*P. 
  }
  else if (isDivisible(qLC, myLC)) { // qLC is divisible by myLC
    LC = div_exact(qLC, myLC);	 //
    if ((LC != 1) && (LC != -1)) { // IMPORTANT: unlike the previous
      // case, we need an additional condition
      // that LC != -1.  THis is because
      // if (LC = -1), then we have qLC and
      // myLC are mutually divisible, and
      // we would be updating R twice!
      R.setCoeff(0, LC); 	   // C = LC, 
      R.setCoeff(R.degree, 1);     // M = X^(myDeg-pDeg)(THIS WAS NOT DONE)
      mulScalar(LC); 	   	   // THIS => THIS * LC

    }
  } else {  			// myLC and qLC are not mutually divisible
    NT g = gcd(qLC, myLC); 	// This ASSUMES gcd is defined in NT !!
    //NT g = 1;  			// In case no gcd is available
    if (g == 1) {
      R.setCoeff(0, qLC);	  	// C = qLC
      R.setCoeff(R.degree, myLC);	 // M = (myLC) * X^{myDeg-pDeg}
      mulScalar(qLC);	 		// forming  C * THIS
      q.mulScalar(myLC);		// forming  M * P
    } else {
      NT qLC1= div_exact(qLC,g);
      NT myLC1= div_exact(myLC,g);
      R.setCoeff(0, qLC1);	  	// C = qLC/g
      R.setCoeff(R.degree, myLC1);	// M = (myLC/g) * X^{myDeg-pDeg}
      mulScalar(qLC1);	 	// forming  C * THIS
      q.mulScalar(myLC1);		// forming  M * P
    }
  }
  (*this) -= q;		// THAT = (C * THIS) - (M * P)

  contract();

  return R;		// Returns R = C + X*M
}// reduceStep

// For internal use only:
// Checks that c*A = B*m + AA 
// 	where A=(*oldthis) and AA=(*newthis)
template <class NT>
Polynomial<NT> Polynomial<NT>::testReduceStep(const Polynomial<NT>& A, 
	const Polynomial<NT>& B) {
std::cout << "+++++++++++++++++++++TEST REDUCE STEP+++++++++++++++++++++\n";
  Polynomial<NT> cA(A);
  Polynomial<NT> AA(A);
  Polynomial<NT> quo;
  quo = AA.reduceStep(B);	        // quo = c + X*m  (m is monomial, c const)
                                // where c*A = B*m + (*newthis)
std::cout << "A = " << A << std::endl;
std::cout << "B = " << B << std::endl;
  cA.mulScalar(quo.coeff[0]);    // A -> c*A
  Polynomial<NT> m(quo);
  m.mulXpower(-1);            // m's value is now m
std::cout << "c = " << quo.coeff[0] << std::endl;
std::cout << "c + xm = " << quo << std::endl;
std::cout << "c*A = " << cA << std::endl;
std::cout << "AA = " << AA << std::endl;
std::cout << "B*m = " << B*m << std::endl;
std::cout << "B*m + AA = " << B*m + AA << std::endl;
  if (cA == (B*m + AA))
	  std::cout << "CORRECT inside testReduceStep" << std::endl;
  else
	  std::cout << "ERROR inside testReduceStep" << std::endl;
std::cout << "+++++++++++++++++END TEST REDUCE STEP+++++++++++++++++++++\n";
  return quo;
}

// PSEUDO-REMAINDER and PSEUDO-QUOTIENT:
// Let A = (*this) and B be the argument polynomial.
// Let Quo be the returned polynomial, 
// and let the final value of (*this) be Rem.
// Also, C is the constant that we maintain.
// We are computing A divided by B.  The relation we guarantee is 
// 		(C * A) = (Quo * B)  + Rem
// where deg(Rem) < deg(B).  So Rem is the Pseudo-Remainder
// and Quo is the Pseudo-Quotient.
// Moreover, C is uniquely determined (we won't spell it out)
// except to note that
//	C divides D = (LC)^{deg(A)-deg(B)+1}
//	where LC is the leading coefficient of B.
// NOTE: 1. Normally, Pseudo-Remainder is defined so that C = D
// 	 So be careful when using our algorithm.
// 	 2. We provide a version of pseudoRemainder which does not
// 	 require C as argument.  [For efficiency, we should provide this
//	 version directly, instead of calling the other version!]

template <class NT>
Polynomial<NT> Polynomial<NT>::pseudoRemainder (
  const Polynomial<NT>& B) {
	NT temp;	// dummy argument to be discarded
	return pseudoRemainder(B, temp);
}//pseudoRemainder

template <class NT>
Polynomial<NT> Polynomial<NT>::pseudoRemainder (
  const Polynomial<NT>& B, NT & C) { 
  contract();         // Let A = (*this).  Contract A.
  Polynomial<NT> tmpB(B);
  tmpB.contract();    // local copy of B
  C = *(new NT(1));  // Initialized to C=1.
  if (B.degree == -1)  {
    std::cout << "ERROR in Polynomial<NT>::pseudoRemainder :\n" <<
    "    -- divide by zero polynomial" << std::endl;
    return Polynomial(0);  // Unit Polynomial (arbitrary!)
  }
  if (B.degree > degree) {
    return Polynomial(); // Zero Polynomial
    // CHECK: 1*THIS = 0*B + THAT,  deg(THAT) < deg(B)
  }

  Polynomial<NT> Quo;  // accumulate the return polynomial, Quo
  Polynomial<NT> tmpQuo;
  while (degree >= B.degree) {  // INVARIANT: C*A = B*Quo + (*this)
    tmpQuo = reduceStep(tmpB);  // Let (*this) be (*oldthis), which
			        // is transformed into (*newthis). Then,
                                //     c*(*oldthis) = B*m + (*newthis)
                                // where tmpQuo = c + X*m
    // Hence,   C*A =   B*Quo +   (*oldthis)      -- the old invariant
    //        c*C*A = c*B*Quo + c*(*oldthis)
    //              = c*B*Quo + (B*m + (*newthis))
    //              = B*(c*Quo + m)  + (*newthis)
    // i.e, to update invariant, we do C->c*C,  Quo --> c*Quo + m.
    C *= tmpQuo.coeff[0];	    // C = c*C
    Quo.mulScalar(tmpQuo.coeff[0]); // Quo -> c*Quo
    tmpQuo.mulXpower(-1);           // tmpQuo is now equal to m
    Quo += tmpQuo;                  // Quo -> Quo + m 
  }

  return Quo;	// Quo is the pseudo-quotient
}//pseudoRemainder

// Returns the negative of the pseudo-remainder
// 	(self-modification)
template <class NT>
Polynomial<NT> & Polynomial<NT>::negPseudoRemainder (
  const Polynomial<NT>& B) {
	NT temp;	// dummy argument to be discarded
	pseudoRemainder(B, temp);
	if (temp < 0) return (*this);
	return negate();
}

template <class NT>
Polynomial<NT> & Polynomial<NT>::operator-() {	// unary minus
  for (int i=0; i<=degree; i++)
    coeff[i] *= -1;
  return *this;
}
template <class NT>
Polynomial<NT> & Polynomial<NT>::power(unsigned int n) {	// self-power
  if (n == 0) {
    degree = 0;
    delete [] coeff;
    coeff = new NT[1];
    coeff[0] = 1;
  } else {
    Polynomial<NT> p = *this;
    for (unsigned int i=1; i<n; i++)
      *this *= p;		// Would a binary power algorithm be better?
  }
  return *this;
}

// evaluation of BigFloat value
//   NOTE: we think of the polynomial as an analytic function in this setting
//      If the coefficients are more general than BigFloats,
//      then we may get inexact outputs, EVEN if the input value f is exact.
//      This is because we would have to convert these coefficients into
//      BigFloats, and this conversion precision is controlled by the
//      global variables defRelPrec and defAbsPrec.
//   
/*
template <class NT>
BigFloat Polynomial<NT>::eval(const BigFloat& f) const {	// evaluation
  if (degree == -1)
    return BigFloat(0);
  if (degree == 0)
    return BigFloat(coeff[0]);
  BigFloat val(0);
  for (int i=degree; i>=0; i--) {
    val *= f;
    val += BigFloat(coeff[i]);	
  }
  return val;
}//eval
*/

// Evaluation Function (generic version)
//
//  NOTE: to call this method, a constructor of the form T(NT) must
//    be available.  Here NT is the type of the coefficients.
//
// ASSERT:
// 		Type T must be >= Type NT
//
// E.g., if NT is BigRat, then T be either BigRat or Expr, but NOT BigFloat.
// 	
// E.g., if NT is BigFloat, it is assumed that the BigFloat has
// 	no error.  Thus, T can be BigFloat, BigRat or Expr in this case.

template <class NT>
template <class T>
T Polynomial<NT>::eval(const T& f) const {	// evaluation
  if (degree == -1)
    return T(0);
  if (degree == 0)
    return T(coeff[0]);
  T val(0);
  for (int i=degree; i>=0; i--) {
    val *= f;
    val += T(coeff[i]);	
  }
  return val;
}//eval


// Evaluation Filter function:
//
// 	ASSERT: NT = BigRat or Expr
//
template <class NT>
BigFloat Polynomial<NT>::evalFilter(const BigFloat& f, bool& validFlag, 
	const extLong& r, const extLong& a) const {	// evaluation
  if (degree == -1)
    return BigFloat(0);
  if (degree == 0)
    return BigFloat(coeff[0], r, a);

  BigFloat fErrFree(f), val(0), valErrFree(0), c;
  fErrFree.makeExact();
  for (int i=degree; i>=0; i--) {
    c = BigFloat(coeff[i], r, a);	
    val *= f; 
    val += c;
    valErrFree *= fErrFree;
    valErrFree += c.makeExact();	
  }
  validFlag = !val.isZeroIn();
  return valErrFree;
}//evalFilter

template <>
CORE_INLINE
BigFloat Polynomial<BigInt>::evalFilter(const BigFloat& f, bool& validFlag, 
	const extLong& r, const extLong& a) const {	// evaluation
  assert(0);
  return BigFloat(0);
}

template <class NT>
BigFloat Polynomial<NT>::evalExact(const Expr& e) const {
  Expr eVal = eval(e);
  eVal.approx();
  return eVal.BigFloatValue();
} 
//============================================================
// Bounds
//============================================================

// Cauchy Upper Bound on Roots
// -- ASSERTION: NT is an integer type
template < class NT >
BigFloat Polynomial<NT>::CauchyUpperBound() const {
  if (zeroP(*this))
    return BigFloat(0);
  NT mx = 0;
  int deg = getTrueDegree();
  for (int i = 0; i < deg; ++i) {
    mx = core_max(mx, abs(coeff[i]));
  }
  Expr e = mx;
  e /= Expr(abs(coeff[deg]));
  e.approx(CORE_INFTY, 2);
  // get an absolute approximate value with error < 1/4
  return (e.BigFloatValue().makeExact() + 2);
}
// Cauchy Lower Bound on Roots
// -- ASSERTION: NT is an integer type
template < class NT >
BigFloat Polynomial<NT>::CauchyLowerBound() const {
  if ((zeroP(*this)) || coeff[0] == 0)
    return BigFloat(0);
  NT mx = 0;
  int deg = getTrueDegree();
  for (int i = 1; i <= deg; ++i) {
    mx = core_max(mx, abs(coeff[i]));
  }
  Expr e = Expr(abs(coeff[0]))/ Expr(abs(coeff[0]) + mx);
  e.approx(2, CORE_INFTY);
  // get an relative approximate value with error < 1/4
  return (e.BigFloatValue().makeExact().div2());
}

// Separation bound for polynomials that may have multiple roots.
// We use the Rump-Schwartz bound.
//
//    ASSERT(the return value is an exact BigFloat and a Lower Bound)
//
template < class NT >
BigFloat Polynomial<NT>::sepBound() const {
  BigInt d;
  BigFloat e;
  int deg = getTrueDegree();

  CORE::power(d, BigInt(deg), ((deg)+4)/2);
  e = CORE::power(height()+1, deg);
  e.makeCeilExact(); // see NOTE below
  return (1/(e*2*d)).makeFloorExact();
        // BUG fix: ``return 1/e*2*d'' was wrong
        // NOTE: the relative error in this division (1/(e*2*d))
        //   is defBFdivRelPrec (a global variable), but
        //   since this is always positive, we are OK.
        //   But to ensure that defBFdivRelPrec is used,
        //   we must make sure that e and d are exact.
        // Finally, using "makeFloorExact()" is OK because
        //   the mantissa minus error (i.e., m-err) will remain positive
        //   as long as the relative error (defBFdivRelPrec) is >1.
}

// height function
template < class NT >
BigFloat Polynomial<NT>::height() const {
  if (zeroP(*this))
    return NT(0);
  int deg = getTrueDegree();
  NT ht = 0;
  for (int i = 0; i< deg; i++)
    if (ht < abs(coeff[i]))
      ht = abs(coeff[i]);
  return BigFloat(ht);
}

// length function
template < class NT >
BigFloat Polynomial<NT>::length() const {
  if (zeroP(*this))
    return NT(0);
  int deg = getTrueDegree();
  NT length = 0;
  for (int i = 0; i< deg; i++)
    length += abs(coeff[i]*coeff[i]);
  return sqrt(BigFloat(length));
}

//============================================================
// differentiation
//============================================================

template <class NT>
Polynomial<NT> & Polynomial<NT>::differentiate() {	// self-differentiation
  if (degree >= 0) {
    NT * c = new NT[degree];
    for (int i=1; i<=degree; i++)
      c[i-1] = coeff[i] * i;
    degree--;
    delete[] coeff;
    coeff = c;
  }
  return *this;
}// differentiation

// multi-differentiate
template <class NT>
Polynomial<NT> & Polynomial<NT>::differentiate(int n) {
  assert(n >= 0);
  for (int i=1; i<=n; i++)
    this->differentiate();
  return *this;
} // multi-differentiate


// ==================================================
// GCD, content, primitive and square-free parts
// ==================================================

/// divisibility predicate for polynomials
// isDivisible(P,Q) returns true iff Q divides P
// To FIX: the predicate name is consistent with Expr.h but not with BigInt.h
template <class NT>
bool isDivisible(Polynomial<NT> p, Polynomial<NT> q) {
  if(zeroP(p))
    return true;
  if(zeroP(q))
    return false;  // We should really return error!