sturm.h

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

      }
    }
  }//isolateRoots(x,y,v)

  // isolateRoots(v)
  ///   isolates all roots and returns them in v
  /**   v is a vector of isolated intervals
   */
  void isolateRoots(BFVecInterval &v) const {
    if (len <= 0) {
       v.clear(); return;
    }
    BigFloat bd = seq[0].CauchyUpperBound();
    // Note: bd is an exact BigFloat (this is important)
    isolateRoots(-bd, bd, v);
  }

  // isolateRoot(i)
  ///   Isolates the i-th smallest root 
  ///         If i<0, isolate the (-i)-th largest root
  ///   Defaults to i=0 (i.e., the smallest positive root a.k.a. main root)
  BFInterval isolateRoot(int i = 0) const {
    if (len <= 0) 
       return BFInterval(1,0);   // ERROR CONDITION
    if (i == 0)
      return mainRoot();
    BigFloat bd = seq[0].CauchyUpperBound();
    return isolateRoot(i, -bd, bd);
  }

  // isolateRoot(i, x, y)
  ///   isolates the i-th smallest root in [x,y]
  /**   If i is negative, then we want the i-th largest root in [x,y]
   *    We assume i is not zero.
   */
  BFInterval isolateRoot(int i, BigFloat x, BigFloat y) const {
    int n = numberOfRoots(x,y);
    if (i < 0) {//then we want the n-i+1 root
      i += n+1;
      if (i <= 0)
        return BFInterval(1,0); // ERROR CONDITION
    }
    if (n < i)
      return BFInterval(1,0);  // ERROR CONDITION INDICATED
    //Now 0< i <= n
    if (n == 1) {
      if ((x>0) || (y<0)) return BFInterval(x, y);
      if (seq[0].coeff[0] == NT(0)) return BFInterval(0,0);
      if (numberOfRoots(0, y)==0) return BFInterval(x,0);
      return BFInterval(0,y);
    }
    BigFloat m = (x+y).div2();
    n = numberOfRoots(x, m);
    if (n >= i)
	    return isolateRoot(i, x, m);
    // Now (n < i) but we have to be careful if m is a root
    if (sign(seq[0].evalExactSign(m)) != 0)   // usual case
      return isolateRoot(i-n, m, y);
    else
      return isolateRoot(i-n+1, m, y);
  }

  // same as isolateRoot(i).
  BFInterval diamond(int i) const {
    return isolateRoot(i);
  }

  // First root above
  BFInterval firstRootAbove(const BigFloat &e) const {
    if (len <= 0)
       return BFInterval(1,0);   // ERROR CONDITION
    return isolateRoot(1, e, seq[0].CauchyUpperBound());
  }

  // Main root (i.e., first root above 0)
  BFInterval mainRoot() const {
    if (len <= 0)
       return BFInterval(1,0);   // ERROR CONDITION
    return isolateRoot(1, 0, seq[0].CauchyUpperBound());
  }

  // First root below
  BFInterval firstRootBelow(const BigFloat &e) const {
    if (len <= 0)
       return BFInterval(1,0);   // ERROR CONDITION
    BigFloat bd = seq[0].CauchyUpperBound(); // bd is exact
    int n = numberOfRoots(-bd, e);
    if (n <= 0)
      return BFInterval(1,0);
    BigFloat bdBF = BigFloat(ceil(bd));
    if (n == 1)
      return BFInterval(-bdBF, e);
    return isolateRoot(n, -bdBF, e);
  }

  // Refine an interval I to absolute precision 2^{-aprec}
  //   THIS USES bisection only!  Use only for debugging (it is too slow)
  //
  BFInterval refine(const BFInterval& I, int aprec) const {
    // assert( There is a unique root in I )
    // We repeat binary search till the following holds
    //      width/2^n <= eps             (eps = 2^(-aprec))
    //   => log(width/eps) <= n
    //   => n = ceil(log(width/eps)) this many steps of binary search
    //   will work.
    // At each step we verify
    //   seq[0].evalExactSign(J.first) * seq[0].evalExactSign(J.second) < 0

    BigFloat width = I.second - I.first;
    if (width <= 0) return I;  // Nothing to do if the
                               //   interval I is exact or inconsistent
    BigFloat eps = BigFloat::exp2(-aprec);   //  eps = 2^{-aprec}
    extLong n =  width.uMSB() + (extLong)aprec;


    BFInterval J = I;           // Return value is the Interval J
    BigFloat midpoint;
    while(n >= 0) {
      midpoint = (J.second + J.first).div2();
      BigFloat m = seq[0].evalExactSign(midpoint);
      if (m == 0) {
        J.first = J.second = midpoint;
        return J;
      }
      if (seq[0].evalExactSign(J.first) * m < 0) {
        J.second = midpoint;
      } else {
        J.first = midpoint;
      }

      n--;
    }

    return J;
  }//End Refine

  // Refine First root above
  BFInterval refinefirstRootAbove(const BigFloat &e, int aprec) const {
    BFInterval I = firstRootAbove(e);
    return refine(I,aprec);
  }

  // Refine First root below
  BFInterval refinefirstRootBelow(const BigFloat &e, int aprec) const {
    BFInterval I = firstRootBelow(e);
    return refine(I,aprec);
  }

  // refineAllRoots(v, aprec)
  //     will modify v so that v is a list of isolating intervals for
  //     the roots of the polynomial in *this.  The size of these intervals
  //     are at most 2^{-aprec}.
  // If v is non-null, we assume it is a list of initial isolating intervals.
  // If v is null, we will first call isolateRoots(v) to set this up.
  void refineAllRoots( BFVecInterval &v, int aprec) {
    BFVecInterval v1;
    BFInterval  J;
    if (v.empty())
      isolateRoots(v);

    for (BFVecInterval::const_iterator it = v.begin();
         it != v.end(); ++it) {        // Iterate through all the intervals
      //refine them to the given precision aprec
      J = refine(BFInterval(it->first, it->second), aprec);
      v1.push_back(std::make_pair(J.first, J.second));
    }
    v.swap(v1);
  }//End of refineAllRoots

  // This is the new version of "refineAllRoots"
  //    	based on Newton iteration
  // It should be used instead of refineAllRoots!
  void newtonRefineAllRoots( BFVecInterval &v, int aprec) {

    BFVecInterval v1;
    BFInterval  J;

    if (v.empty())
      isolateRoots(v);
    for (BFVecInterval::iterator it = v.begin();
         it != v.end(); ++it) {        // Iterate through all the intervals
      //refine them to the given precision aprec
      J = newtonRefine(*it, aprec);

      if (NEWTON_DIV_BY_ZERO) {
        J.first = 1;
        J.second = 0;   // indicating divide by zero
      }
      v1.push_back(std::make_pair(J.first, J.second));
    }
    v.swap(v1);
  }//End of newtonRefineAllRoots

  /** val = newtonIterN(n, bf, del, err, fuMSB, ffuMSB)
   * 
   *    val is the root after n iterations of Newton
   *       starting from initial value of bf and is exact.
   *    fuMSB and ffuMSB are precision parameters for the approximating
   *		the coefficients of the underlyinbg polynomial, f(x).
   *    	THEY are used ONLY if the coefficients of the polynomial
   *		comes from a field (in particular, Expr or BigRat).
   *		We initially approximate the coefficients of f(x) to fuMSB 
   *		relative bits, and f'(x) to ffuMSB relative bits.
   *		The returned values of fuMSB and ffuMSB are the final
   *		precision used by the polynomial evaluation algorithm.
   *    Return by reference, "del" (difference between returned val and value
   *       in the previous Newton iteration)
   *
   *    Also, "err" is returned by reference and bounds the error in "del".
   *
   *    IMPORTANT: we assume that when x is an exact BigFloat,
   *    then Polynomial<NT>::eval(x) will be exact!
   *    But current implementation of eval() requires NT <= BigFloat.
   * ****************************************************/    

  BigFloat newtonIterN(long n, const BigFloat& bf, BigFloat& del,
	unsigned long & err, extLong& fuMSB, extLong& ffuMSB) {
    if (len <= 0) return bf;   // Nothing to do!  User must
                               // check this possibility!
    BigFloat val = bf;  
    // val.makeExact();    // val is exact

    // newton iteration
    for (int i=0; i<n; i++) {
      ////////////////////////////////////////////////////
      // Filtered Eval
      ////////////////////////////////////////////////////
      BigFloat ff = seq[1].evalExactSign(val, 3*ffuMSB); //3 is a slight hack
      ffuMSB = ff.uMSB();
      //ff is guaranteed to have the correct sign as the exact evaluation.
      ////////////////////////////////////////////////////

      if (ff == 0) {
        NEWTON_DIV_BY_ZERO = true;
        del = 0;
        core_error("Zero divisor in Newton Iteration",
                __FILE__, __LINE__, false);
        return 0;
      }

      ////////////////////////////////////////////////////
      // Filtered Eval
      ////////////////////////////////////////////////////
      BigFloat f= seq[0].evalExactSign(val, 3*fuMSB); //3 is a slight hack
      fuMSB = f.uMSB();
      ////////////////////////////////////////////////////

      if (f == 0) {
        NEWTON_DIV_BY_ZERO = false;
        del = 0;    // Indicates that we have reached the exact root
		    //    This is because eval(val) is exact!!!
        return val; // val is the exact root, before the last iteration
      }
      del = f/ff; // But the accuracy of "f/ff" must be controllable
		    // by the caller...
      err = del.err();
      del.makeExact(); // makeExact() is necessary
      val -= del;
      // val.makeExact();  // -- unnecessary...
    }
    return val;
  }//newtonIterN

  //Another version of newtonIterN which does not return the error 
  //and passing the uMSB as arguments; it is easier for the user to call
  //this.
  BigFloat newtonIterN(long n, const BigFloat& bf, BigFloat& del){
    unsigned long err;
    extLong fuMSB=0, ffuMSB=0;
    return newtonIterN(n, bf, del, err, fuMSB, ffuMSB);
  }

  // v = newtonIterE(prec, bf, del, fuMSB, ffuMSB)
  //
  //    return the value v which is obtained by Newton iteration
  //    until del.uMSB < -prec, starting from initial value of bf.
  //    Returned value is an exact BigFloat.
  //    We guarantee at least one Newton step (so del is defined).
  //
  //	   The parameters fuMSB and ffuMSB are precision parameters for
  //	   evaluating coefficients of f(x) and f'(x), used similarly
  //	   as described above for newtonIterN(....)
  //
  //    Return by reference "del" (difference between returned val and value
  //       in the previous Newton iteration).  This "del" is an upper bound
  //       on the last (f/f')-value in Newton iteration.
  //
  //    IN particular, if v is in the Newton zone of a root z^*, then z^* is
  //       guaranteed to lie inside [v-del, v+del].
  //
  //    Note that this is dangerous unless you know that bf is already
  //       in the Newton zone.  So we use the global N_STOP_ITER to
  //       prevent infinite loop.

  BigFloat newtonIterE(int prec, const BigFloat& bf, BigFloat& del, 
	extLong& fuMSB, extLong& ffuMSB) {
    // usually, prec is positive
    int count = N_STOP_ITER; // upper bound on number of iterations
    int stepsize = 1;
    BigFloat val = bf;
    unsigned long err = 0;

    do {
      val = newtonIterN(stepsize, val, del, err, fuMSB, ffuMSB);
      count -= stepsize;
      stepsize++; // heuristic
    } while ((del != 0) && ((del.uMSB() >= -prec) && (count >0))) ;

    if (count == 0) core_error("newtonIterE: reached count=0",
		    	__FILE__, __LINE__, true);
    del = BigFloat(core_abs(del.m()), err, del.exp() );
    del.makeCeilExact();
    return val;
  }

  //Another version of newtonIterE which avoids passing the uMSB's.
  BigFloat newtonIterE(int prec, const BigFloat& bf, BigFloat& del){
    extLong fuMSB=0, ffuMSB=0;
    return newtonIterE(prec, bf, del, fuMSB, ffuMSB);
  }
  // A Smale bound which is an \'a posteriori condition. Applying 
  // Newton iteration to any point z satisfying this condition we are 
  // sure to converge to the nearest root in a certain interval of z.
  // The condition is for all k >= 2,
  //    | \frac{p^(k)(z)}{k!p'(z)} |^{1\(k-1)} < 1/8 * |\frac{p'(z)}{p(z)}|
  // Note: code below has been streamlined (Chee)
  /*
    bool smaleBound(const Polynomial<NT> * p, BigFloat z){
    int deg = p[0].getTrueDegree();
    BigFloat max, temp, temp1, temp2;
    temp2 = p[1].eval(z);
    temp = core_abs(temp2/p[0].eval(z))/8;
    BigInt fact_k = 2;
    for(int k = 2; k <= deg; k++){
      temp1 = core_abs(p[k].eval(z)/(fact_k*temp2)); 
      if(k-1 == 2)
	temp1 = sqrt(temp1);
      else
	temp1 = root(temp1, k-1);
      if(temp1 >= temp) return false; 
    }
    return true;
    }
   */

  //An easily computable Smale's point estimate for Newton as compared to the
  //one above. The criterion is
  //
  // ||f||_{\infty} * \frac{|f(z)|}{|f'(z)|^2} 
  //                * \frac{\phi'(|z|)^2}{\phi(|z|)}  < 0.03
  // where
  //           \phi(r) = \sum_{i=0}{m}r^i,
  //           m = deg(f)
  //
  //It is given as Theorem B in [Smale86].
  //Reference:- Chapter 8 in Complexity and Real Computation, by
  //            Blum, Cucker, Shub and Smale
  //
  //For our implementation we calculate an upper bound on
  //the second fraction in the inequality above.  For r>0,
  //
  //    \phi'(r)^2     m^2 (r^m + 1)^2
  //     ---------  <  -------------------          (1)
  //    \phi(r)        (r-1) (r^{m+1} - 1)
  //
  // Alternatively, we have
  // 
  //    \phi'(r)^2     (mr^{m+1} + 1)^2
  //     ---------  <  -------------------          (2)
  //    \phi(r)        (r-1)^3 (r^{m+1} - 1)
  //
  // The first bound is better when r > 1.
  // The second bound is better when r << 1.
  // Both bounds (1) and (2) assumes r is not equal to 1.
  // When r=1, the exact value is
  //
  //    \phi'(r)^2     m^2 (m + 1)
  //     ---------  =  -----------                  (3)
  //    \phi(r)            4
  // 
  // REMARK: smaleBoundTest(z) actually computes an upper bound
  // 	on alpha(f,z), and compares it to 0.02 (then our theory
  // 	says that z is a robust approximate zero).