sturm.h

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

  //
  bool smaleBoundTest(const BigFloat& z){
    assert(z.isExact());   // the bound only makes sense for exact z

#ifdef CORE_DEBUG
    std::cout <<"Computing Smale's bound = " <<  std::endl;
#endif

    if(seq[0].evalExactSign(z) == 0)// Reached the exact root.
      return true;

    BigFloat fprime = core_abs(seq[1].evalExactSign(z));
    fprime.makeFloorExact();
    if (fprime == 0) return false;  // z is a critical value!
    BigFloat temp =        // evalExactSign(z) may have error.
      core_abs(seq[0].evalExactSign(z));
    temp = (temp.makeCeilExact()/power(fprime, 2)).makeCeilExact();
    temp = temp*seq[0].height();  // remains exact
    //Thus, temp >=  ||f||_{\infty} |\frac{f(z)}{f'(z)^2}|

    int m = seq[0].getTrueDegree();    
    BigFloat x = core_abs(z);
    if (x==1)   // special case, using (3)
	    return (temp * BigFloat(m*m*(m+1)).div2().div2() < 0.02);

    BigFloat temp1;
    if (x>1) { // use formula (1)
      temp1 = power(m* (power(x, m)+1), 2);          // m^2*(x^m + 1)^2
      temp1 /= ((x - 1)*(power(x, m+1) - 1));        // formula (1)
    } else {  // use formula (2)
      temp1 = power(m*(power(x, m+1) +1), 2);        // (m*x^{m+1} + 1)^2
      temp1 /= (power(x - 1,3)*(power(x, m+1) -1));  // formula (2)
    }

#ifdef CORE_DEBUG
    std::cout <<"Value returned by Smale bound = " << temp * temp1.makeCeilExact() << std::endl;
#endif

    if(temp * temp1.makeCeilExact() < 0.03)          // make temp1 exact!
      return true;
    else
      return false;
  }//smaleBoundTest


  // yapsBound(p)
  // 	returns a bound on size of isolating interval of polynomial p
  // 	which is guaranteed to be in the Newton Zone.
  //    N.B. p MUST be square-free
  //
  //   Reference: Theorem 6.37, p.184 of Yap's book
  //   	   [Fundamental Problems of Algorithmic Algebra]

  BigFloat yapsBound(const Polynomial<NT> & p) const {
    int deg = p.getTrueDegree();
    return  1/(1 + pow(BigFloat(deg), 3*deg+9)
               *pow(BigFloat(2+p.height()),6*deg));
  }

  //newtonRefine(J, a) 
  //
  //    ASSERT(J is an isolating interval for some root x^*)
  //
  //    ASSERT(J.first and J.second are exact BigFloats)
  //
  //    Otherwise, the boundaries of the interval are not well defined.
  //    We will return a refined interval with exact endpoints,
  //    still called J, containing x^* and
  //
  // 			|J| < 2^{-a}.
  //
  // 	TO DO: write a version of newtonRefine(J, a, sign) where
  // 	sign=J.first.sign(), since you may already know the sign
  // 	of J.first.  This will skip the preliminary stuff in the
  // 	current version.
  //
  BFInterval newtonRefine(BFInterval &J, int aprec) {

#ifdef CORE_DEBUG_NEWTON
std::cout << "In newtonRefine, input J=" << J.first
	<< ", " << J.second << " precision = " << aprec << std::endl;
#endif

    if (len <= 0) return J;   // Nothing to do!  User must
                               // check this possibility!
      

    if((J.second - J.first).uMSB() < -aprec){
      return (J);
    }
    int xSign, leftSign, rightSign;

    leftSign = sign(seq[0].evalExactSign(J.first));
    if (leftSign == 0) {
      J.second = J.first;
      return J;
    }

    rightSign = sign(seq[0].evalExactSign(J.second));
    if (rightSign == 0) {
      J.first = J.second;
      return J;
    }

    assert( leftSign * rightSign < 0 );

    //N is number of times Newton is called without checking
    // whether the result is still in the interval or not
    #define NO_STEPS 2
    // REMARK: NO_STEPS=1 is incorrect, as it may lead to
    //      linear convergence (it is somewhat similar to Dekker-Brent's
    //      idea of guaranteeing that bisection does not
    //	    destroy the superlinear convergence of Newton.
    int N = NO_STEPS;

    BigFloat x, del, olddel, temp;
    unsigned long err;
    BigFloat yap = yapsBound(seq[0]);

    BigFloat old_width = J.second - J.first;
    x = (J.second + J.first).div2();

    // initial estimate for the evaluation of filter to floating point precision
    extLong fuMSB=54, ffuMSB=54;

    //MAIN WHILE LOOP. We ensure that J always contains the root

    while ( !smaleBoundTest(x) && 
	    (J.second - J.first) > yap &&
	   (J.second - J.first).uMSB() >= -aprec) {
      x = newtonIterN(N, x, del, err, fuMSB, ffuMSB);
      if ((del == 0)&&(NEWTON_DIV_BY_ZERO == false)) {  // reached exact root!
        J.first = J.second = x;
        return J;
      }

      BigFloat left(x), right(x);
      if (del>0) {
      	left -= del; right += del;
      } else {
      	left += del; right -= del;
      }

      // update interval
      if ((left > J.first)&&(left <J.second)) {
	  int lSign = sign(seq[0].evalExactSign(left));
          if (lSign == leftSign)  // leftSign=sign of J.first
            J.first = left;
	  else if (lSign == 0) {
            J.first = J.second = left;
            return J;
          } else {
	    J.second = left;
          }	
      }
      if ((right < J.second)&&(right >J.first)) {
	  int rSign = sign(seq[0].evalExactSign(right));
          if (rSign == rightSign)
            J.second = right;
	  else if (rSign == 0) {
            J.first = J.second = right;
            return J;
          } else {
            J.first = right;
          }
      }
      BigFloat width = J.second - J.first;

      //left and right are exact, since x is exact.
      if (width*2 <= old_width && !NEWTON_DIV_BY_ZERO) {
                                  // we can get a better root:

	// No, it is not necessary to update x to
	// the midpoint of the new interval J.
	// REASON: basically, it is hard to be smarter than Newton's method!
	// Newton might bring x very close to one endpoint, but it can be
	// because the root is near there!  In any case,
	// by setting x to the center of J, you only gain at most
	// one bit of accuracy, but you stand to loose an
	// arbitrary amount of bits of accuracy if you are unlucky!
	// So I will comment out the next line.  --Chee (Aug 9, 2004).
	// 
	// x = (J.second + J.first).div2();
	if (J.first > x || J.second < x)
	  x = (J.second + J.first).div2();

	old_width = width; // update width

        N ++;      // be more confident or aggressive
	           //  (perhaps we should double N)
		   //
      } else {// Either NEWTON_DIV_BY_ZERO=true
	      // Or width has not decreased sufficiently
	x = (J.second + J.first).div2();//Reset x to midpoint since it was the
	                                //value from a failed Newton step
	xSign = sign(seq[0].evalExactSign(x));
	if (xSign == rightSign) {
	  J.second = x;
	} else if (xSign == leftSign) {
	  J.first = x;
	} else { // xSign must be 0
	  J.first = J.second = x; return J;
	}
	x = (J.second + J.first).div2();

	old_width = old_width.div2(); // update width
	
	// reduce value of N:
	N = core_max(N-1, NO_STEPS);   // N must be at least NO_STEPS
      }
    }//MAIN WHILE LOOP

    if((J.second - J.first).uMSB() >= -aprec){ // The interval J
                    //still hasn't reached the required precision.
                    //But we know the current value of x (call it x_0)
		    //is in the strong Newton basin of the
		    //root x^* (because it passes Smale's bound)
      //////////////////////////////////////////////////////////////////
      //Both x_0 and the root x^* are in the interval J.
      //Let NB(x^*) be the strong Newton basin of x^*.  By definition,
      //x_0 is in NB(x^*) means that:
      //
      //    x_0 is in NB(x^*) iff |x_n-x^*| \le 2^{-2^{n}+1} |x_0-x^*|
      //    
      // where x_n is the n-th iterate of Newton.  
      //    
      //  LEMMA 1: if x_0  \in NB(x^*) then 
      //               |x_0 - x^*| <= 2|del|      (*)
      //  and
      //               |x_1 - x^*| <= |del|       (**)
      //
      //  where del = -f(x_0)/f'(x_0) and x_1 = x_0 + del
      //Proof:
      //Since x_0 is in the strong Newton basin, we have
      //         |x_1-x^*| <= |x_0-x^*|/2.        (***)
      //The bound (*) is equivalent to
      //         |x_0-x^*|/2 <= |del|.
      //This is equivalent to
      //         |x_0-x^*| - |del| <= |x_0-x^*|/2,
      //which follows from
      //         |x_0-x^* + del| <= |x_0-x^*|/2,
      //which is equivalent to (***).  
      //The bound (**) follows from (*) and (***).
      //QED
      //
      //  COMMENT: the above derivation refers to the exact values,
      //  but what we compute is X_1 where X_1 is an approximation to
      //  x_1.  However, let us write X_1 = x_0 - DEL, where DEL is
      //  an approximation to del.  
      //
      //  LEMMA 2:  If |DEL| >= |del|,
      //  then (**) holds with X_1 and DEL in place of x_1 and del.
      //  
      //  NOTE: We implemented this DEL in newtonIterE.   

#ifdef CORE_DEBUG
      std::cout << "Inside Newton Refine: Refining Part " << std::endl;

      if((J.second - J.first) > yap)
	std::cout << "Smales Bound satisfied " << std::endl;
      else
	std::cout << "Chees Bound satisfied " << std::endl;
#endif
      xSign = sign(seq[0].evalExactSign(x));
      if(xSign == 0){
	J.first = J.second = x; 
	return J; // missing before!
      }

      //int k = clLg((-(J.second - J.first).lMSB() + aprec).asLong());
      x = newtonIterE(aprec, x, del, fuMSB, ffuMSB);
      xSign = sign(seq[0].evalExactSign(x));

      if(xSign == leftSign){//Root is greater than x
	J.first = x;
	J.second = x + del;  // justified by Lemma 2 above
      }else if(xSign == rightSign){//Root is less than x
	J.first = x - del;   // justified by Lemma 2 above
	J.second = x ;
      }else{//x is the root
	J.first = J.second = x;
      }
    }



#ifdef CORE_DEBUG
    std::cout << " Returning from Newton Refine: J.first = " << J.first
	      << " J.second = " << J.second << " aprec = " << aprec
	      << " Sign at the interval endpoints = " 
	      << sign(seq[0].evalExactSign(J.first))
	      << " : " << sign(seq[0].evalExactSign(J.second)) << " Err at starting = " 
	      << J.first.err() << " Err at end = " << J.second.err() << std::endl;
#endif

    assert( (seq[0].evalExactSign(J.first) * seq[0].evalExactSign(J.second) <= 0) );

#ifdef CORE_DEBUG_NEWTON
    if (seq[0].evalExactSign(J.first) * seq[0].evalExactSign(J.second) > 0)
      std::cout <<" ERROR! Root is not in the Interval " << std::endl;
    if(J.second - J.first >  BigFloat(1).exp2(-aprec))
      std::cout << "ERROR! Newton Refine failed to achieve the desired precision" << std::endl;
#endif

      return(J);
 }//End of newton refine

};// Sturm class

// ==================================================
// Static initialization
// ==================================================
template <class NT>
int Sturm<NT>:: N_STOP_ITER = 10000;   // stop IterE after this many loops
// Reset this as needed

// ==================================================
// Helper Functions 
// ==================================================

// isZeroIn(I):
//          returns true iff 0 is in the closed interval I
CORE_INLINE bool isZeroIn(BFInterval I) {
	return ((I.first <= 0.0) && (I.second >= 0.0));
}

/////////////////////////////////////////////////////////////////
//  DIAGNOSTIC TOOLS
/////////////////////////////////////////////////////////////////
// Polynomial tester:   P is polynomial to be tested
//          prec is the bit precision for root isolation
//          n is the number of roots predicted

template<class NT>
CORE_INLINE void testSturm(const Polynomial<NT>&P, int prec, int n = -1) {
  Sturm<NT> SS (P);
  BFVecInterval v;
  SS.refineAllRoots(v, prec);
  std::cout << "   Number of roots is " << v.size() <<std::endl;
  if ((n >= 0) & (v.size() == (unsigned)n))
    std::cout << " (CORRECT!)" << std::endl;
  else
    std::cout << " (ERROR!) " << std::endl;
  int i = 0;
  for (BFVecInterval::const_iterator it = v.begin();
       it != v.end(); ++it) {
    std::cout << ++i << "th Root is in ["
    << it->first << " ; " << it->second << "]" << std::endl;
  }
}// testSturm

// testNewtonSturm( Poly, aprec, n)
//   will run the Newton-Sturm refinement to isolate the roots of Poly
//         until absolute precision aprec.
//   n is the predicated number of roots
//      (will print an error message if n is wrong)
template<class NT>
CORE_INLINE void testNewtonSturm(const Polynomial<NT>&P, int prec, int n = -1) {
  Sturm<NT> SS (P);
  BFVecInterval v;
  SS.newtonRefineAllRoots(v, prec);
  std::cout << "   Number of roots is " << v.size();
  if ((n >= 0) & (v.size() == (unsigned)n))
    std::cout << " (CORRECT!)" << std::endl;
  else
    std::cout << " (ERROR!) " << std::endl;

  int i = 0;
  for (BFVecInterval::iterator it = v.begin();
       it != v.end(); ++it) {
    std::cout << ++i << "th Root is in ["
    << it->first << " ; " << it->second << "]" << std::endl;
    if(it->second - it->first <= (1/power(BigFloat(2), prec)))
      std::cout << " (CORRECT!) Precision attained" << std::endl;
    else
      std::cout << " (ERROR!) Precision not attained" << std::endl;
  }
}// testNewtonSturm

CORE_END_NAMESPACE

#endif