poly.tcc

很多二维 三维几何计算算法 C++ 类库

// 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, always returns the exact value)
///
///  This evaluation function is easy to use, but may not be efficient
///  when you have BigRat or Expr values.
///
/// User must be aware that the return type of eval is Max of Types NT and T.
///
/// E.g., If NT is BigRat, and T is Expr then Max(NT,T)=Expr. 
/// 	
/// REMARK: If NT is BigFloat, it is assumed that the BigFloat is error-free.  

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


/// Approximate Evaluation of Polynomials
/// 	the coefficients of the polynomial are approximated to some
///	specified composite precision (r,a).
/// @param  f evaluation point 
/// @param  r relative precision to which the coefficients are evaluated
/// @param  a absolute precision to which the coefficients are evaluated
/// @return a BigFloat with error containing value of the polynomial.
///     If zero is in this BigFloat interval, then we don't know the sign.
//
// 	ASSERT: NT = BigRat or Expr
//
template <class NT>
BigFloat Polynomial<NT>::evalApprox(const BigFloat& f, 
	const extLong& r, const extLong& a) const {	// evaluation
  if (degree == -1)
    return BigFloat(0);
  if (degree == 0)
    return BigFloat(coeff[0], r, a);

  BigFloat val(0), c;
  for (int i=degree; i>=0; i--) {
    c = BigFloat(coeff[i], r, a);	
    val *= f; 
    val += c;
  }
  return val;
}//evalApprox

// This BigInt version of evalApprox should never be called...
template <>
CORE_INLINE
BigFloat Polynomial<BigInt>::evalApprox(const BigFloat& /*f*/,
	const extLong& /*r*/, const extLong& /*a*/) const {	// evaluation
  assert(0);
  return BigFloat(0);
}



/**
 * Evaluation at a BigFloat value
 * using "filter" only when NT is BigRat or Expr.
 * Using filter means we approximate the polynomials
 * coefficients using BigFloats.  If this does not give us
 * the correct sign, we will resort to an "exact" evaluation
 * using Expr.
 *
 * If NT <= BigFloat, we just call eval().
 *
   We use the following heuristic estimates of precision for coefficients:

      r = 1 + lg(|P|_\infty) + lg(d+1)  		if f <= 1
      r = 1 + lg(|P|_\infty) + lg(d+1) + d*lg|f| 	if f > 1
      
   if the filter fails, then we use Expr to do evaluation.

   This function is mainly called by Newton iteration (which
   has some estimate for oldMSB from previous iteration).

   @param p polynomial to be evaluated
   @param val the evaluation point
   @param oldMSB an rough estimate of the lg|p(val)|
   @return bigFloat interval contain p(val), with the correct sign

 ***************************************************/
template <class NT>
BigFloat Polynomial<NT>::evalExactSign(const BigFloat& val,
	 const extLong& oldMSB) const {
    assert(val.isExact());
    if (getTrueDegree() == -1)
      return BigFloat(0);
  
    extLong r;
    r = 1 + BigFloat(height()).uMSB() + clLg(long(getTrueDegree()+1));
    if (val > 1)
      r += getTrueDegree() * val.uMSB();
    r += core_max(extLong(0), -oldMSB);
  
    if (hasExactDivision<NT>::check()) { // currently, only to detect NT=Expr and NT=BigRat
        BigFloat rVal = evalApprox(val, r);
        if (rVal.isZeroIn()) {
	  Expr eVal = eval(Expr(val));	// eval gives exact value
	  eVal.approx(54,CORE_INFTY);  // if value is 0, we get exact 0
	  return eVal.BigFloatValue();
	} else 
          return rVal;
    } else
	return BigFloat(eval(val));

   //return 0; // unreachable
  }//evalExactSign
  

//============================================================
// 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);
}

//============================================================
// An iterative version of computing Cauchy Bound from Erich Kaltofen.
// See the writeup under collab/sep/.
//============================================================
template < class NT >
BigInt Polynomial<NT>::CauchyBound() const {
  int deg = getTrueDegree();
  BigInt B(1);
  BigFloat lhs(0), rhs(1);
  while (true) {
    /* compute \sum_{i=0}{deg-1}{|a_i|B^i} */
    lhs = 0;
    for (int i=deg-1; i>=0; i--) {
      lhs *= B;
      lhs += abs(coeff[i]);
    }
    lhs /= abs(coeff[deg]);
    lhs.makeFloorExact();
    /* compute B^{deg} */
    if (rhs <= lhs) {
      B <<= 1;
      rhs *= (BigInt(1)<<deg);
    } else
      break;
  }
  return B;
}

//====================================================================
//Another iterative bound which is at least as good as the above bound
//by Erich Kaltofen.
//====================================================================
template < class NT >
BigInt Polynomial<NT>::UpperBound() const {
  int deg = getTrueDegree();

  BigInt B(1);
  BigFloat lhsPos(0), lhsNeg(0), rhs(1);
  while (true) {
    /* compute \sum_{i=0}{deg-1}{|a_i|B^i} */
    lhsPos = lhsNeg = 0;
    for (int i=deg-1; i>=0; i--) {
      if (coeff[i]>0) {
      	lhsPos = lhsPos * B + coeff[i];
      	lhsNeg = lhsNeg * B;
      } else {
      	lhsNeg = lhsNeg * B - coeff[i];
      	lhsPos = lhsPos * B;
      } 
    }
    lhsNeg /= abs(coeff[deg]);
    lhsPos /= abs(coeff[deg]);
    lhsPos.makeCeilExact();
    lhsNeg.makeCeilExact();

    /* compute B^{deg} */
    if (rhs <= max(lhsPos,lhsNeg)) {
      B <<= 1;
      rhs *= (BigInt(1)<<deg);
    } else
      break;
  }
  return B;
}

// 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
/// @return a BigFloat with error
template < class NT >
BigFloat Polynomial<NT>::height() const {
  if (zeroP(*this))
    return BigFloat(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
/// @return a BigFloat with error
template < class NT >
BigFloat Polynomial<NT>::length() const {
  if (zeroP(*this))
    return BigFloat(0);
  int deg = getTrueDegree();
  NT length = 0;
  for (int i = 0; i< deg; i++)
    length += abs(coeff[i]*coeff[i]);
  return sqrt(BigFloat(length));
}