sturm.h

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

/****************************************************************************
 * Core Library Version 1.7, August 2004
 * Copyright (c) 1995-2004 Exact Computation Project
 * All rights reserved.
 *
 * This file is part of CORE (http://cs.nyu.edu/exact/core/); you may
 * redistribute it under the terms of the Q Public License version 1.0.
 * See the file LICENSE.QPL distributed with CORE.
 *
 * Licensees holding a valid commercial license may use this file in
 * accordance with the commercial license agreement provided with the
 * software.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 *  File: Sturm.h
 * 
 *  Description: 
 *  The templated class Sturm implements Sturm sequences.
 *  Basic capabilities include:
 *     counting number of roots in an interval, 
 *     isolating all roots in an interval
 *     isolating the i-th largest (or smallest) root in interval
 *  It is based on the Polynomial class.
 * 
 *   BigFloat intervals are used for this (new) version.
 *   It is very important that the BigFloats used in these intervals
 *   have no error at the beginning, and this is maintained
 *   by refinement.  Note that if x, y are error-free BigFloats,
 *   then (x+y)/2 may not be error-free (in current implementaion.
 *   We have to call a special "exact divide by 2" method,
 *   (x+y).div2() for this purpose.
 *
 *   CONVENTION: an interval defined by a pair of BigFloats x, y
 *   has this interpretation:
 *       (1) if x>y,  it represents an invalid interval.
 *       (2) if x=y,  it represents a unique point x.
 *       (3) if x<y,  it represents the open interval (x,y).
 *           In this case, we always may sure that x, y are not zeros.
 * 
 *   TODO LIST and Potential Bugs:
 *   (1) Split an isolating interval to give definite sign (done)
 *   (2) Should have a test for square-free polynomials (done)
 * 
 *  Author:  Chee Yap and Sylvain Pion, Vikram Sharma
 *  Date:    July 20, 2002
 *
 * WWW URL: http://cs.nyu.edu/exact/
 * Email: [e-mail]
 *
 * $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.3-branch/Core/include/CGAL/CORE/poly/Sturm.h $
 * $Id: Sturm.h 37060 2007-03-13 18:10:39Z reichel $
 ***************************************************************************/


#ifndef CORE_STURM_H
#define CORE_STURM_H

#include "CGAL/CORE/BigFloat.h"
#include "CGAL/CORE/Expr.h"
#include "CGAL/CORE/poly/Poly.h"

CORE_BEGIN_NAMESPACE

// ==================================================
// Sturm Class
// ==================================================

template < class NT >
class Sturm {
public:
  int len;      // len is 1 less than the number of non-zero entries in array seq.
  		//     I.e., len + 1 = length of the Sturm Sequence
                // N.B. When len = -1 or len = 0 are special,
                //     the array seq is not used!
                //     Hence, one must test these special cases
  Polynomial<NT> * seq;      // array of polynomials of length "len+1"
  Polynomial<NT> g;//GCD of input polynomial P and it's derivative P'
  NT cont;//Content of the square-free part of input polynomial P
  //Thus P = g * cont * seq[0]
  static int N_STOP_ITER;    // Stop IterE after this many iterations. This
  // is initialized below, outside the Newton class
  bool NEWTON_DIV_BY_ZERO;   // This is set to true when there is divide by
  // zero in Newton iteration (at critical value)
  // User is responsible to check this and to reset.
  typedef Polynomial<NT> PolyNT;

  // ===============================================================
  // CONSTRUCTORS
  // ===============================================================
  // Null Constructor
  Sturm() : len(0), NEWTON_DIV_BY_ZERO(false) {}

  // Constructor from a polynomial
  Sturm(Polynomial<NT> pp) : NEWTON_DIV_BY_ZERO(false) {
    len = pp.getTrueDegree();
    if (len <= 0) return; // hence, seq is not defined in these cases
    seq = new Polynomial<NT> [len+1];
    seq[0] = pp;
    g = seq[0].sqFreePart();
    cont = content(seq[0]);
    seq[0].primPart();
    seq[1] = differentiate(seq[0]);
    int i;
    for (i=2; i <= len; i++) {
      seq[i] = seq[i-2];
      seq[i].negPseudoRemainder(seq[i-1]);
      if (zeroP(seq[i])){
	len = i-1;//Since len is one less than the number of non-zero entries.
        break;
      }
      seq[i].primPart(); // Primitive part is important to speed
      // up large polynomials! However, for first 2 polymials,
      // we MUST NOT take primitive part, because we
      // want to use them in Newton Iteration
    }
  }

  // Chee: 7/31/04
  // 	We need BigFloat version of Sturm(Polynomial<NT>pp) because
  // 	of curve verticalIntersection() ... .  We also introduce
  // 	various support methods in BigFloat.h (exact_div, gcd, etc).
  // Constructor from a BigFloat polynomial 
  //	Need the fake argument to avoid compiler overloading errors
  Sturm(Polynomial<BigFloat> pp, bool fake) : NEWTON_DIV_BY_ZERO(false) {
    len = pp.getTrueDegree();
    if (len <= 0) return; // hence, seq is not defined in these cases
    seq = new Polynomial<NT> [len+1];
    seq[0] = pp;
    g = seq[0].sqFreePart();
    cont = content(seq[0]);
    seq[0].primPart();
    seq[1] = differentiate(seq[0]);
    int i;
    for (i=2; i <= len; i++) {
      seq[i] = seq[i-2];
      seq[i].negPseudoRemainder(seq[i-1]);
      if (zeroP(seq[i])){
	len = i-1;//Since len is one less than the number of non-zero entries.
	//len = i;
        break;
      }
      seq[i].primPart(); // Primitive part is important to speed
      // up large polynomials! However, for first 2 polymials,
      // we DO NOT take primitive part, because we
      // want to use them in Newton Iteration
    }
  }

  // Constructor from an array of NT's
  //   -- this code is identical to constructing from a polynomial...
  Sturm(int n, NT * c) : NEWTON_DIV_BY_ZERO(false) {
    Polynomial<NT> pp(n, c); // create the polynomial pp first and call the
    (*this) = Sturm<NT>(pp);//constructor from a polynomial
  }

  // copy constructor
  Sturm(const Sturm&s) : len(s.len), NEWTON_DIV_BY_ZERO(s.NEWTON_DIV_BY_ZERO) {
    if (len <= 0) return;
    seq = new Polynomial<NT> [len+1];
    for (int i=0; i<=len; i++)
      seq[i] = s.seq[i];
  }

  // assignment operator
  const Sturm& operator=(const Sturm& o) {
    if (this == &o)
      return *this;
    if (len > 0)
      delete[] seq;
    NEWTON_DIV_BY_ZERO = o.NEWTON_DIV_BY_ZERO;
    len = o.len;
    if (len > 0) {
      seq = new Polynomial<NT>[len+1];
      for (int i=0; i<=len; i++)
        seq[i] = o.seq[i];
    }
    return *this;
  }

  // destructor
  ~Sturm() {
    if (len != 0)
      delete[] seq;
  }

  // METHODS

  // dump functions
  void dump(std::string msg) const {
    std::cerr << msg << std::endl;
    if (len <= 0) std::cerr << " len = " << len << std::endl;
    else
       for (int i=0; i<=len; i++)
         std::cerr << " seq[" << i << "] = " << seq[i] << std::endl;
  }
  void dump() const {
    dump("");
  }

  // signVariations(x, sx)
  //   where sx = sign of evaluating seq[0] at x
  //   PRE-CONDITION: sx != 0  and len > 0
  int signVariations(const BigFloat & x, int sx) const {
    assert((sx != 0) && (len >0));
    int cnt = 0;
    int last_sign = sx;
    for (int i=1; i<=len; i++) {// Chee (4/29/04): Bug fix,
        // should start iteration at i=1, not i=0.  Potential error
        // if seq[0].eval(x)=0 (though not in our usage).
      int sgn = sign(seq[i].evalExactSign(x));
      if (sgn*last_sign < 0) {
        cnt++;
        last_sign *= -1;
      }
    }
    return cnt;
  }

  // signVariations(x)
  //   --the first polynomial eval is not yet done
  //   --special return value of -1, indicating x is root!
  int signVariations(const BigFloat & x) const {
    if (len <= 0) return len; 
    int signx = sign(seq[0].evalExactSign(x));
    if (signx == 0)
      return (-1);    // THIS indicates that x is a root...
    		      // REMARK: in our usage, this case does not arise
    return signVariations(x, signx);
  }//signVariations(x)

  // signVariation at +Infinity
  int signVariationsAtPosInfty() const {
    if (len <= 0) return len;
    int cnt = 0;
    int last_sign = sign(seq[0].coeff[seq[0].getTrueDegree()]);
    assert(last_sign != 0);
    for (int i=1; i<=len; i++) {
      int sgn = sign(seq[i].coeff[seq[i].getTrueDegree()]);
      if (sgn*last_sign < 0)
        cnt++;
      if (sgn != 0)
        last_sign = sgn;
    }
    return cnt;
  }

  // signVariation at -Infinity
  int signVariationsAtNegInfty() const {
    if (len <= 0) return len;
    int cnt = 0;
    int last_sign = sign(seq[0].coeff[seq[0].getTrueDegree()]);
    if (seq[0].getTrueDegree() % 2 != 0)
      last_sign *= -1;
    assert(last_sign != 0);
    for (int i=1; i<=len; i++) {
      int parity = (seq[i].getTrueDegree() % 2 == 0) ? 1 : -1;
      int sgn = parity * sign(seq[i].coeff[seq[i].getTrueDegree()]);
      if (sgn*last_sign < 0)
        cnt++;
      if (sgn != 0)
        last_sign = sgn;
    }
    return cnt;
  }

  // numberOfRoots(x,y):
  //   COUNT NUMBER OF ROOTS in the close interval [x,y]
  //   IMPORTANT: Must get it right even if x, y are roots
  //   Assert("x and y are exact")
  //       [If the user is unsure of this assertion, do
  //        "x.makeExact(); y.makeExact()" before calling].
  ///////////////////////////////////////////
  int numberOfRoots(const BigFloat &x, const BigFloat &y) const {
    assert(x <= y);   // we allow x=y
    if (len <= 0) return len;  // return of -1 means infinity of roots!
    int signx = sign(seq[0].evalExactSign(x));
    if (x == y) return ((signx == 0) ? 1 : 0);
    int signy = sign(seq[0].evalExactSign(y));
    // easy case: THIS SHOULD BE THE OVERWHELMING MAJORITY

    if (signx != 0 && signy != 0)
      return (signVariations(x, signx) - signVariations(y, signy));
    // harder case: THIS SHOULD BE VERY INFREQUENT
    BigFloat sep = (seq[0].sepBound()).div2();
    BigFloat newx, newy;
    if (signx == 0)
      newx = x - sep;
    else
      newx = x;
    if (signy == 0)
      newy = y + sep;
    else
      newy = y;
    return (signVariations(newx, sign(seq[0].evalExactSign(newx)))
            - signVariations(newy, sign(seq[0].evalExactSign(newy))) );
  }//numberOfRoots

  // numberOfRoots():
  //   Counts the number of real roots of a polynomial
  ///////////////////////////////////////////
  int numberOfRoots() const {
    if (len <= 0) return len;  // return of -1 means infinity of roots!
    //    BigFloat bd = seq[0].CauchyUpperBound();
    //    return numberOfRoots(-bd, bd);
    return signVariationsAtNegInfty() - signVariationsAtPosInfty();
  }

  // numberOfRoots above or equal to x:
  //   Default value x=0 (i.e., number of positive roots)
  //   assert(len >= 0)
  ///////////////////////////////////////////
  int numberOfRootsAbove(const BigFloat &x = 0) const {
    if (len <= 0) return len;  // return of -1 means infinity of roots!
    int signx = sign(seq[0].evalExactSign(x));
    if (signx != 0)
      return signVariations(x, signx) - signVariationsAtPosInfty();
    BigFloat newx = x - (seq[0].sepBound()).div2();
    return signVariations(newx, sign(seq[0].evalExactSign(newx)))
           - signVariationsAtPosInfty();
  }

  // numberOfRoots below or equal to x:
  //   Default value x=0 (i.e., number of negative roots)
  //   assert(len >= 0)
  ///////////////////////////////////////////
  int numberOfRootsBelow(const BigFloat &x = 0) const {
    if (len <= 0) return len;  // return of -1 means infinity of roots!
    int signx = sign(seq[0].evalExactSign(x));
    if (signx != 0)
      return signVariationsAtNegInfty() - signVariations(x, signx);
    BigFloat newx = x + (seq[0].sepBound()).div2();
    return signVariationsAtNegInfty()
           - signVariations(newx, sign(seq[0].evalExactSign(newx)));
  }


  /// isolateRoots(x, y, v)
  ///             Assertion(x, y are exact BigFloats)
  ///   isolates all the roots in [x,y] and returns them in v.
  /**   v is a list of intervals
   *    [x,y] is the initial interval to be isolated
   *
   *    Properties we guarantee in the return values:
   *
   *    (0) All the intervals have exact BigFloats as endpoints
   *    (1) If 0 is a root, the corresponding isolating interval will be
   *        exact, i.e., we return [0,0].
   *    (2) If an interval is [0,x], it contains a positive root
   *    (3) If an interval is [y,0], it contains a negative root
   */
  void isolateRoots(const BigFloat &x, const BigFloat &y,
                    BFVecInterval &v) const {
    assert(x<=y);

    int n = numberOfRoots(x,y);
    if (n == 0) return;
    if (n == 1) {
      if ((x > 0) || (y < 0)) // usual case: 0 is not in interval
        v.push_back(std::make_pair(x, y));
      else { // if 0 is inside our interval (this extra
	     // service is not strictly necessary!)
        if (seq[0].coeff[0] == 0)
          v.push_back(std::make_pair(BigFloat(0), BigFloat(0)));
        else if (numberOfRoots(0,y) == 0)
          v.push_back(std::make_pair(x, BigFloat(0)));
        else
          v.push_back(std::make_pair(BigFloat(0), y));
      }
    } else { // n > 1
      BigFloat mid = (x+y).div2(); // So mid is exact.
      if (sign(seq[0].evalExactSign(mid)) != 0)  { // usual case: mid is non-root
      	isolateRoots(x, mid, v);
      	isolateRoots(mid, y, v); 
      } else { // special case: mid is a root
	BigFloat tmpEps = (seq[0].sepBound()).div2();  // this is exact!
	if(mid-tmpEps > x )//Since otherwise there are no roots in (x,mid)
	  isolateRoots(x, (mid-tmpEps).makeCeilExact(), v);
	v.push_back(std::make_pair(mid, mid));
	if(mid+tmpEps < y)//Since otherwise there are no roots in (mid,y)
	  isolateRoots((mid+tmpEps).makeFloorExact(), y, v);