exprrep.h

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

  void reduceToBigRat(const BigRat&);
  /// reduce current node
  void reduceTo(const ExprRep*);
  /// reduce current node to zero
  void reduceToZero();
  /// return operator string
  virtual const std::string op() const {
    return "UNKNOWN";
  }
  //@}

  /// \name Degree Bound Functions
  //@{
  /// compute "d_e" based on # of sqrts
  extLong degreeBound();
  /// count actually computes the degree bound of current node.
  virtual extLong count() = 0;
  /// reset the flag "visited"
  virtual void clearFlag() = 0;
  //@}
#ifdef CORE_DEBUG
  virtual unsigned long dagSize() = 0;
  virtual void fullClearFlag() = 0;
#endif
};//ExprRep

/// \class ConstRep
/// \brief constant node
class ConstRep : public ExprRep {
public:
  /// \name Constructors and Destructor
  //@{
  /// default constructor
  ConstRep() {}
  /// destructor
  virtual ~ConstRep() {}
  //@}

  /// \name Debug Functions
  //@{
  /// print debug information in list mode
  void debugList(int level, int depthLimit) const;
  /// print debug information in tree mode
  void debugTree(int level, int indent, int depthLimit) const;
  //@}
protected:
  /// initialize nodeInfo
  virtual void initNodeInfo();
  /// return operator in string
  const std::string op() const {
    return "C";
  }
  /// count returns the degree of current node
  //extLong count() { return d_e(); }
  extLong count();
  /// clear visited flag
  void clearFlag() {
    visited() = false;
  }
#ifdef CORE_DEBUG
  unsigned long dagSize();
  void fullClearFlag();
#endif
};

/// \class ConstDoubleRep
/// \brief constant node
class ConstDoubleRep : public ConstRep {
public:
  /// \name Constructors and Destructor
  //@{
  /// default constructor
  ConstDoubleRep() {}
  /// constructor for all \c double type
  ConstDoubleRep(double d) {
    ffVal = d;
  }
  /// destructor
  ~ConstDoubleRep() {}
  //@}
  CORE_MEMORY(ConstDoubleRep)
protected:
  /// compute sign and MSB
  void computeExactFlags();
  /// compute approximation value
  void computeApproxValue(const extLong&, const extLong&);
};

/// \class ConstRealRep
/// \brief constant node
class ConstRealRep : public ConstRep {
public:
  /// \name Constructors and Destructor
  //@{
  /// default constructor
  ConstRealRep() : value(CORE_REAL_ZERO) { }
  /// constructor for all \c Real type
  ConstRealRep(const Real &);
  /// destructor
  ~ConstRealRep() {}
  //@}
  CORE_MEMORY(ConstRealRep)
private:
  Real value; ///< internal representation of node
protected:
  /// compute sign and MSB
  void computeExactFlags();
  /// compute approximation value
  void computeApproxValue(const extLong&, const extLong&);
};

/// \class Constant Polynomial Node
/// \brief template class where NT is supposed to be some number type
template <class NT>
class ConstPolyRep : public ConstRep {
public:
  /// \name Constructors and Destructor
  //@{
  /// default constructor
  ConstPolyRep() { }

  /// constructor for Polynomial
  ConstPolyRep(const Polynomial<NT>& p, int n) : ss(p) {
    // isolate roots using Sturm Sequences
    I = ss.isolateRoot(n);
    // check whether n-th root exists
    if (I.first == 1 && I.second == 0) {
      core_error("CORE ERROR! root index out of bound",
		      __FILE__, __LINE__, true);
      abort();
    }
    // test if the root isolated in I is 0:
    if ((I.first == 0)&&(I.second == 0))
      ffVal = 0;
    else
      ffVal = computeFilteredValue();	// silly to use a filter here!
    					// since sign is known.
  }

  /// constructor for Polynomial
  ConstPolyRep(const Polynomial<NT>& p, const BFInterval& II)
	  : ss(p), I(II) {
    BFVecInterval v;
    ss.isolateRoots(I.first, I.second, v);
    I = v.front();
    if (v.size() != 1) {
      core_error("CORE ERROR! non-isolating interval",
		      __FILE__, __LINE__, true);
      abort();
    }
    ffVal = computeFilteredValue(); // Chee: this line seems unnecessary
  }

  /// destructor
  ~ConstPolyRep() {}
  //@}
  CORE_MEMORY(ConstPolyRep)

private:
  Sturm<NT> ss; ///< internal Sturm sequences
  BFInterval I; ///< current interval contains the real value
  		// IMPORTANT: I.first and I.second are exact BigFloats
  filteredFp computeFilteredValue() {
    // refine initial interval to absolute error of 2^(lMSB(k)-54) where
    // 	  k is a lower bound on the root (use Cauchy Lower Bound).
    //    Hence, the precision we pass to refine should be 54-lMSB(k).

    // refine with newton (new method)
    // ss.seq[0] could be zero!!
    // I=ss.newtonRefine(I,
    //            54-(ss.seq[0].CauchyLowerBound()).lMSB().asLong());
    extLong lbd = ss.seq[0].CauchyLowerBound().lMSB();

    if (lbd.isTiny())
      I = ss.newtonRefine(I, 54);
    else
      I = ss.newtonRefine(I, 54-lbd.asLong()); // is this necessary?

    //return I.first.doubleValue(); // NOTE:  This is not quite right!
    // It should be "centralized" to set
    // the error bit correctly.
    // E.g., otherwise, radical(4,2) will print wrongly.
    if ((I.first == 0) && (I.second == 0))	// Checkfor zero value
      return filteredFp(0);
    BigFloat x = centerize(I.first, I.second);
    double val = x.doubleValue();
    double max = core_max(core_abs(I.first), core_abs(I.second)).doubleValue();
    int ind = 1;
    /*
    long ee = x.exp()*CHUNK_BIT;
    unsigned long err = ee > 0 ? (x.err() << ee) : (x.err() >> (-ee));
    double max = core_abs(val) + err;
    int ind = longValue((BigInt(x.err()) << 53) / (x.m() + x.err())); 
    */
    return filteredFp(val, max, ind); // Aug 8, 2004, Comment from Chee:
       // I think we should get rid of filters here!  Given the interval I,
       // we either know the sign (I.first >=0) or (I.second <=0)
       // or we don't.  We don't need to compute all the index stuff.
       // In fact, you have lost the sign in the above computation...
       // ALSO, why bother to use filter?
  }//computeFilteredValue

protected:
  void initNodeInfo() {
    nodeInfo = new NodeInfo();
    d_e() = ss.seq[0].getTrueDegree(); // return degree of the polynomial
  }
  /// compute sign and MSB
  void computeExactFlags() {

    if ((I.first == 0) && (I.second == 0)) {
      reduceToZero();
      return;
    } else if (I.second > 0) {
      uMSB() = I.second.uMSB();
      lMSB() = I.first.lMSB();
      sign() = 1;
    } else { // we know that I.first < 0
      lMSB() = I.second.lMSB();
      uMSB() = I.first.uMSB();
      sign() = -1;
    }
    // length() = 1+ ss.seq[0].length().uMSB();
    measure() = 1+ ss.seq[0].length().uMSB();	// since measure<= length

    // compute u25, l25, v2p, v2m, v5p, v5m
    v2p() = v2m() = v5p() = v5m() = 0;
    u25() = 1+ss.seq[0].CauchyUpperBound().uMSB();
    l25() = ceilLg(ss.seq[0].getLeadCoeff());  // assumed coeff is integer!!
    		// ceilLg(BigInt) and ceilLg(Expr) are defined. But if
    		// NT=int, ceilLg(int) is ambiguous!  Added ceilLg(int)
		// under BigInt.h

    // compute high, low, lc, tc
    high() = u25();
    low() = - (ss.seq[0].CauchyLowerBound().lMSB()); // note the use of negative
    lc() = l25();
    tc() = ceilLg(ss.seq[0].getTailCoeff());

    // no rational reduction
    if (rationalReduceFlag)
      ratFlag() = -1;

    flagsComputed() = true;
    appValue()=centerize(I.first, I.second);// set an initial value for appValue
  }
  /// compute approximation value
  void computeApproxValue(const extLong& relPrec, const extLong& absPrec) {
    extLong pr = -lMSB() + relPrec;
    extLong p = pr < absPrec ? pr : absPrec;

    // bisection sturm (old method)
    //I = ss.refine(I, p.asLong()+1);

    // refine with newton (new method)
    I = ss.newtonRefine(I, p.asLong()+1);
    appValue() = centerize(I.first, I.second);
  }
};
/// \class UnaryOpRep
/// \brief unary operator node
class UnaryOpRep : public ExprRep {
public:
  /// \name Constructors and Destructor
  //@{
  /// constructor
  UnaryOpRep(ExprRep* c) : child(c)  {
    child->incRef();
  }
  /// destructor
  virtual ~UnaryOpRep() {
    child->decRef();
  }
  //@}

  /// \name Debug Functions
  //@{
  /// print debug information in list mode
  void debugList(int level, int depthLimit) const;
  /// print debug information in tree mode
  void debugTree(int level, int indent, int depthLimit) const;
  //@}
protected:
  ExprRep* child; ///< pointer to its child node
  /// initialize nodeInfo
  virtual void initNodeInfo();
  /// clear visited flag
  void clearFlag();
#ifdef CORE_DEBUG
  unsigned long dagSize();
  void fullClearFlag();
#endif
};

/// \class NegRep
/// \brief unary minus operator node
class NegRep : public UnaryOpRep {
public:
  /// \name Constructors and Destructor
  //@{
  /// constructor
  NegRep(ExprRep* c) : UnaryOpRep(c) {
    ffVal = - child->ffVal;
  }
  /// destructor
  ~NegRep() {}
  //@}

  CORE_MEMORY(NegRep)
protected:
  /// compute sign and MSB
  void computeExactFlags();
  /// compute approximation value
  void computeApproxValue(const extLong&, const extLong&);
  /// return operator in string
  const std::string op() const {
    return "Neg";
  }
  /// count computes the degree of current node, i.e., d_e().
  /** This is now a misnomer, but historically accurate.
   */
  extLong count();
};

/// \class SqrtRep
/// \brief squartroot operator node
class SqrtRep : public UnaryOpRep {
public:
  /// \name Constructors and Destructor
  //@{
  /// constructor
  SqrtRep(ExprRep* c) : UnaryOpRep(c) {
    ffVal = (child->ffVal).sqrt();
  }
  /// destructor
  ~SqrtRep() {}
  //@}

  CORE_MEMORY(SqrtRep)
protected:
  /// compute sign and MSB
  void computeExactFlags();
  /// compute approximation value
  void computeApproxValue(const extLong&, const extLong&);
  /// return operator in string
  const std::string op() const {
    return "Sqrt";
  }
  /// count computes the degree of current node, i.e., d_e().
  /** This is now a misnomer, but historically accurate.
   */
  extLong count();
};

/// \class BinOpRep
/// \brief binary operator node
class BinOpRep : public ExprRep {
public:
  /// \name Constructors and Destructor
  //@{
  /// constructor
  BinOpRep(ExprRep* f, ExprRep* s) : first(f), second(s) {
    first->incRef();
    second->incRef();
  }
  /// destructor
  virtual ~BinOpRep() {
    first->decRef();
    second->decRef();
  }
  //@}

  /// \name Debug Functions
  //@{
  /// print debug information in list mode
  void debugList(int level, int depthLimit) const;
  /// print debug information in tree mode
  void debugTree(int level, int indent, int depthLimit) const;
  //@}
protected:
  ExprRep* first;  ///< first operand
  ExprRep* second; ///< second operand

  /// initialize nodeInfo
  virtual void initNodeInfo();
  /// clear visited flags
  void clearFlag();
  /// count computes the degree of current node, i.e., d_e().
  /** This is now a misnomer, but historically accurate.
   */
  extLong count();
#ifdef CORE_DEBUG
  unsigned long dagSize();
  void fullClearFlag();
#endif
};

/// \struct Add
/// \brief "functor" class used as parameter to AddSubRep<>
struct Add {
  /// name
  static const char* name;

  /// unary operator
  template <class T>
  const T& operator()(const T& t) const {
    return t;
  }

  /// binary operator
  template <class T>
  T operator()(const T& a, const T& b) const {
    return a+b;
  }
};

/// \struct Sub
/// \brief "functor" class used as parameter to AddSubRep<>
struct Sub {
  /// name
  static const char* name;

  /// unary operator
  template <class T>
  T operator()(const T& t) const {
    return -t;
  }

  /// binary operator
  template <class T>
  T operator()(const T& a, const T& b) const {
    return a-b;
  }
};

/// \class AddSubRep
/// \brief template class where operator is supposed to be Add or Sub
template <class Operator>