SPAR/AIModule/BWTA/vendors/CGAL/CGAL/CORE/ExprRep.h

Go to the documentation of this file.

/****************************************************************************
 * 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: ExprRep.h
 * Synopsis: Internal Representation of Expr.
 * 
 * Written by 
 *       Koji Ouchi <ouchi@simulation.nyu.edu>
 *       Chee Yap <yap@cs.nyu.edu>
 *       Igor Pechtchanski <pechtcha@cs.nyu.edu>
 *       Vijay Karamcheti <vijayk@cs.nyu.edu>
 *       Chen Li <chenli@cs.nyu.edu>
 *       Zilin Du <zilin@cs.nyu.edu>
 *       Sylvain Pion <pion@cs.nyu.edu> 
 *       Vikram Sharma<sharma@cs.nyu.edu>
 *
 * WWW URL: http://cs.nyu.edu/exact/
 * Email: exact@cs.nyu.edu
 *
 * $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.5-branch/Core/include/CGAL/CORE/ExprRep.h $
 * $Id: ExprRep.h 41739 2008-01-20 23:09:58Z spion $
 ***************************************************************************/

#ifndef _CORE_EXPRREP_H_
#define _CORE_EXPRREP_H_

#include <CGAL/CORE/Real.h>
#include <CGAL/CORE/Filter.h>
#include <CGAL/CORE/poly/Sturm.h>

CORE_BEGIN_NAMESPACE

#ifdef CORE_DEBUG_BOUND
// These counters are incremented each time each bound is recognized as equal
// to the best one in computeBound().
extern unsigned int BFMSS_counter;
extern unsigned int Measure_counter;
// extern unsigned int Cauchy_counter;
extern unsigned int LiYap_counter;
// These counters are incremented each time each bound is recognized as equal
// to the best one in computeBound(), and it's strictly the best.
extern unsigned int BFMSS_only_counter;
extern unsigned int Measure_only_counter;
// extern unsigned int Cauchy_only_counter;
extern unsigned int LiYap_only_counter;
// This counter is incremented each time the precision needed matches the
// root bound.
extern unsigned int rootBoundHitCounter;
#endif

const extLong EXTLONG_BIG = (1L << 30);
const extLong EXTLONG_SMALL = -(1L << 30);

const double log_5 = log(double(5))/log(double(2));

// Returns the ceil of log_2(5^a).
inline extLong ceilLg5(const extLong & a) {
#if defined (_MSC_VER)
  return (int) ::ceil(log_5 * a.toLong());
#else
  return (int) std::ceil(log_5 * a.toLong());
#endif
}

struct NodeInfo {
  Real     appValue;            
  bool     appComputed;         
  bool     flagsComputed;       
  extLong  knownPrecision;      

#ifdef CORE_DEBUG
  extLong relPrecision;
  extLong absPrecision;
  unsigned long numNodes;
#endif


  extLong d_e;

  bool visited;         
  int sign;             

  extLong  uMSB; 
  extLong  lMSB; 

  // For the degree-length method mentioned in Chee's book.
  /* the degree of defining polynomial P(X) obtained from Resultant calculus
   * (deprecated now) */
  // extLong degree;

  // extLong length; ///< length is really lg(|| P(X) ||)
  extLong measure; 

  // For our new bound.

  extLong high;

  extLong low;
  extLong lc;
  extLong tc;

  // For the 2-ary BFMSS bound.
  extLong v2p, v2m;
  // For the 5-ary BFMSS bound.
  extLong v5p, v5m;

  extLong u25;

  extLong l25;

  int ratFlag;          
  BigRat* ratValue;     

  NodeInfo();
};//NodeInfo struct

//  forward reference
// class Expr;

//  Members: private: int refCount,
//            public:  NodeInfo* nodeInfo,
//                     filteredFp ffVal.
class ExprRep {
public:


  ExprRep();
  virtual ~ExprRep() {
    if (nodeInfo != NULL) // This check is only for optimization.
      delete nodeInfo;
  }



  void incRef() {
    ++refCount;
  }
  void decRef() {
    if (--refCount == 0)
      delete this;
  }
  int getRefCount() const {
    return refCount;
  }
  int isUnique() const {
    return refCount == 1;
  }



  const Real & getAppValue(const extLong& relPrec = defRelPrec,
                           const extLong& absPrec = defAbsPrec);
  int getSign();
  int getExactSign();

  const Real& appValue() const {
    return nodeInfo->appValue;
  }
  Real& appValue() {
    return nodeInfo->appValue;
  }

  const bool& appComputed() const {
    return nodeInfo->appComputed;
  }
  bool& appComputed() {
    return nodeInfo->appComputed;
  }

  const bool& flagsComputed() const {
    return nodeInfo->flagsComputed;
  }
  bool& flagsComputed() {
    return nodeInfo->flagsComputed;
  }

  const extLong& knownPrecision() const {
    return nodeInfo->knownPrecision;
  }
  extLong& knownPrecision() {
    return nodeInfo->knownPrecision;
  }

#ifdef CORE_DEBUG
  const extLong& relPrecision() const {
    return nodeInfo->relPrecision;
  }
  extLong& relPrecision() {
    return nodeInfo->relPrecision;
  }

  const extLong& absPrecision() const {
    return nodeInfo->absPrecision;
  }
  extLong& absPrecision() {
    return nodeInfo->absPrecision;
  }

  const unsigned long& numNodes() const {
    return nodeInfo->numNodes;
  }
  unsigned long& numNodes() {
    return nodeInfo->numNodes;
  }
#endif

  const extLong& d_e() const {
    return nodeInfo->d_e;
  }
  extLong& d_e() {
    return nodeInfo->d_e;
  }

  const bool& visited() const {
    return nodeInfo->visited;
  }
  bool& visited() {
    return nodeInfo->visited;
  }

  const int& sign() const {
    return nodeInfo->sign;
  }
  int& sign() {
    return nodeInfo->sign;
  }

  const extLong& uMSB() const {
    return nodeInfo->uMSB;
  }
  extLong& uMSB() {
    return nodeInfo->uMSB;
  }

  const extLong& lMSB() const {
    return nodeInfo->lMSB;
  }
  extLong& lMSB() {
    return nodeInfo->lMSB;
  }

  //  const extLong& length() const { return nodeInfo->length; }
  //  extLong& length() { return nodeInfo->length; }

  const extLong& measure() const {
    return nodeInfo->measure;
  }
  extLong& measure() {
    return nodeInfo->measure;
  }

  const extLong& high() const {
    return nodeInfo->high;
  }
  extLong& high() {
    return nodeInfo->high;
  }

  const extLong& low() const {
    return nodeInfo->low;
  }
  extLong& low() {
    return nodeInfo->low;
  }

  const extLong& lc() const {
    return nodeInfo->lc;
  }
  extLong& lc() {
    return nodeInfo->lc;
  }

  const extLong& tc() const {
    return nodeInfo->tc;
  }
  extLong& tc() {
    return nodeInfo->tc;
  }

  const extLong& v2p() const {
    return nodeInfo->v2p;
  }
  extLong& v2p() {
    return nodeInfo->v2p;
  }

  const extLong& v2m() const {
    return nodeInfo->v2m;
  }
  extLong& v2m() {
    return nodeInfo->v2m;
  }

  extLong v2() const {
    return v2p()-v2m();
  }

  const extLong& v5p() const {
    return nodeInfo->v5p;
  }
  extLong& v5p() {
    return nodeInfo->v5p;
  }

  const extLong& v5m() const {
    return nodeInfo->v5m;
  }
  extLong& v5m() {
    return nodeInfo->v5m;
  }

  extLong v5() const {
    return v5p()-v5m();
  }

  const extLong& u25() const {
    return nodeInfo->u25;
  }
  extLong& u25() {
    return nodeInfo->u25;
  }

  const extLong& l25() const {
    return nodeInfo->l25;
  }
  extLong& l25() {
    return nodeInfo->l25;
  }

  const int& ratFlag() const {
    return nodeInfo->ratFlag;
  }
  int& ratFlag() {
    return nodeInfo->ratFlag;
  }

  const BigRat* ratValue() const {
    return nodeInfo->ratValue;
  }
  BigRat*& ratValue() {
    return nodeInfo->ratValue;
  }

  BigInt BigIntValue();
  BigRat BigRatValue();
  BigFloat BigFloatValue();
  // toString() Joaquin Grech 31/5/2003
  std::string toString(long prec=defOutputDigits, bool sci=false) {
    return (getAppValue(defRelPrec, defAbsPrec)).toString(prec,sci);
  }



  const std::string dump(int = OPERATOR_VALUE) const;
  virtual void debugList(int level, int depthLimit) const = 0;
  virtual void debugTree(int level, int indent, int depthLimit) const = 0;


  friend std::ostream& operator<<(std::ostream&, ExprRep&);

private:
  int refCount;    // reference count

public:
  enum {OPERATOR_ONLY, VALUE_ONLY, OPERATOR_VALUE, FULL_DUMP};

  NodeInfo* nodeInfo;   
  filteredFp ffVal;     



  virtual void initNodeInfo() = 0;
  virtual void computeExactFlags() = 0;
  extLong computeBound();
  void approx(const extLong& relPrec, const extLong& absPrec);
  virtual void computeApproxValue(const extLong&, const extLong&) = 0;
  bool withinKnownPrecision(const extLong&, const extLong&);



  void reduceToBigRat(const BigRat&);
  void reduceTo(const ExprRep*);
  void reduceToZero();
  virtual const std::string op() const {
    return "UNKNOWN";
  }



  extLong degreeBound();
  virtual extLong count() = 0;
  virtual void clearFlag() = 0;
#ifdef CORE_DEBUG
  virtual unsigned long dagSize() = 0;
  virtual void fullClearFlag() = 0;
#endif
};//ExprRep

class ConstRep : public ExprRep {
public:


  ConstRep() {}
  virtual ~ConstRep() {}



  void debugList(int level, int depthLimit) const;
  void debugTree(int level, int indent, int depthLimit) const;
protected:
  virtual void initNodeInfo();
  const std::string op() const {
    return "C";
  }
  //extLong count() { return d_e(); }
  extLong count();
  void clearFlag() {
    visited() = false;
  }
#ifdef CORE_DEBUG
  unsigned long dagSize();
  void fullClearFlag();
#endif
};

class ConstDoubleRep : public ConstRep {
public:


  ConstDoubleRep() {}
  ConstDoubleRep(double d) {
    ffVal = d;
  }
  ~ConstDoubleRep() {}
  CORE_MEMORY(ConstDoubleRep)
protected:
  void computeExactFlags();
  void computeApproxValue(const extLong&, const extLong&);
};

class ConstRealRep : public ConstRep {
public:


  ConstRealRep() : value(CORE_REAL_ZERO) { }
  ConstRealRep(const Real &);
  ~ConstRealRep() {}
  CORE_MEMORY(ConstRealRep)
private:
  Real value; 
protected:
  void computeExactFlags();
  void computeApproxValue(const extLong&, const extLong&);
};

template <class NT>
class ConstPolyRep : public ConstRep {
public:


  ConstPolyRep() { }

  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.
  }

  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
  }

  ~ConstPolyRep() {}
  CORE_MEMORY(ConstPolyRep)

private:
  Sturm<NT> ss; 
  BFInterval I; 
                // 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
  }
  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
  }
  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 : public ExprRep {
public:


  UnaryOpRep(ExprRep* c) : child(c)  {
    child->incRef();
  }
  virtual ~UnaryOpRep() {
    child->decRef();
  }



  void debugList(int level, int depthLimit) const;
  void debugTree(int level, int indent, int depthLimit) const;
protected:
  ExprRep* child; 

  virtual void initNodeInfo();
  void clearFlag();
#ifdef CORE_DEBUG
  unsigned long dagSize();
  void fullClearFlag();
#endif
};

class NegRep : public UnaryOpRep {
public:


  NegRep(ExprRep* c) : UnaryOpRep(c) {
    ffVal = - child->ffVal;
  }
  ~NegRep() {}

  CORE_MEMORY(NegRep)
protected:
  void computeExactFlags();
  void computeApproxValue(const extLong&, const extLong&);
  const std::string op() const {
    return "Neg";
  }

  extLong count();
};

class SqrtRep : public UnaryOpRep {
public:


  SqrtRep(ExprRep* c) : UnaryOpRep(c) {
    ffVal = (child->ffVal).sqrt();
  }
  ~SqrtRep() {}

  CORE_MEMORY(SqrtRep)
protected:
  void computeExactFlags();
  void computeApproxValue(const extLong&, const extLong&);
  const std::string op() const {
    return "Sqrt";
  }

  extLong count();
};

class BinOpRep : public ExprRep {
public:


  BinOpRep(ExprRep* f, ExprRep* s) : first(f), second(s) {
    first->incRef();
    second->incRef();
  }
  virtual ~BinOpRep() {
    first->decRef();
    second->decRef();
  }



  void debugList(int level, int depthLimit) const;
  void debugTree(int level, int indent, int depthLimit) const;
protected:
  ExprRep* first;  
  ExprRep* second; 

  virtual void initNodeInfo();
  void clearFlag();

  extLong count();
#ifdef CORE_DEBUG
  unsigned long dagSize();
  void fullClearFlag();
#endif
};

struct Add {
  static const char* name;

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

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

struct Sub {
  static const char* name;

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

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

template <class Operator>
class AddSubRep : public BinOpRep {
public:


  AddSubRep(ExprRep* f, ExprRep* s) : BinOpRep(f, s) {
    ffVal = Op(first->ffVal, second->ffVal);
  }
  ~AddSubRep() {}

  CORE_MEMORY(AddSubRep)
protected:
  void computeExactFlags();
  void computeApproxValue(const extLong&, const extLong&);
  const std::string op() const {
    return Operator::name;
  }
private:
  static Operator Op;
};//AddSubRep class

template <class Operator>
Operator AddSubRep<Operator>::Op;

template <class Operator>
void AddSubRep<Operator>::computeExactFlags() {
  if (!first->flagsComputed())
    first->computeExactFlags();
  if (!second->flagsComputed())
    second->computeExactFlags();

  int sf = first->sign();
  int ss = second->sign();

  if ((sf == 0) && (ss == 0)) { // the node is zero
    reduceToZero();
    return;
  } else if (sf == 0) { // first operand is zero
    reduceTo(second);
    sign() = Op(ss);
    appValue() = Op(appValue());
    if (rationalReduceFlag && ratFlag() > 0)
      *(ratValue()) = Op(*(ratValue()));
    return;
  } else if (ss == 0) { // second operand is zero
    reduceTo(first);
    return;
  }
  // rational node
  if (rationalReduceFlag) {
    if (first->ratFlag() > 0 && second->ratFlag() > 0) {
      BigRat val=Op(*(first->ratValue()), *(second->ratValue()));
      reduceToBigRat(val);
      ratFlag() = first->ratFlag() + second->ratFlag();
      return;
    } else
      ratFlag() = -1;
  }

  // neither operand is zero
  extLong df    = first->d_e();
  extLong ds    = second->d_e();
  // extLong md    = df < ds ? df : ds;
  // extLong l1    = first->length();
  // extLong l2    = second->length();
  extLong m1    = first->measure();
  extLong m2    = second->measure();

  // length() = df * l2 + ds * l1 + d_e() + md;
  measure() = m1 * ds + m2 * df + d_e();

  // BFMSS[2,5] bound.
  v2p() = core_min(first->v2p() + second->v2m(),
                  first->v2m() + second->v2p());
  v2m() = first->v2m() + second->v2m();
  v5p() = core_min(first->v5p() + second->v5m(),
                  first->v5m() + second->v5p());
  v5m() = first->v5m() + second->v5m();

  if (v2p().isInfty() || v5p().isInfty())
    u25() = CORE_INFTY;
  else
    u25() = EXTLONG_ONE + core_max(first->v2p() + second->v2m()
            - v2p() + ceilLg5(first->v5p() + second->v5m() - v5p())
            + first->u25() + second->l25(),
                                   first->v2m() + second->v2p() - v2p()
            + ceilLg5(first->v5m() + second->v5p() - v5p())
            + first->l25() + second->u25());
  l25() = first->l25() + second->l25();

  lc() = ds * first->lc() + df * second->lc();
  tc() = measure();

  high() = core_max(first->high(),second->high())+EXTLONG_ONE;
  // The following is a subset of the minimization in computeBound().
  low() = core_min(measure(), (d_e()-EXTLONG_ONE)*high() + lc());

  extLong lf = first->lMSB();
  extLong ls = second->lMSB();
  extLong uf = first->uMSB();
  extLong us = second->uMSB();

  extLong l  = core_max(lf, ls);
  extLong u  = core_max(uf, us);

#ifdef CORE_TRACE
  std::cout << "INSIDE Add/sub Rep: " << std::endl;
#endif

  if (Op(sf, ss) != 0) {     // can't possibly cancel out
#ifdef CORE_TRACE
    std::cout << "Add/sub Rep:  Op(sf, ss) non-zero" << std::endl;
#endif

    uMSB() = u + EXTLONG_ONE;
    lMSB() = l;            // lMSB = core_min(lf, ls)+1 better
    sign() = sf;
  } else {               // might cancel out
#ifdef CORE_TRACE
    std::cout << "Add/sub Rep:  Op(sf, ss) zero" << std::endl;
#endif

    uMSB() = u + EXTLONG_ONE;
    uMSB() = u;
    if (lf >= us + EXTLONG_TWO) {// one is at least 1 order of magnitude larger
#ifdef CORE_TRACE
      std::cout << "Add/sub Rep:  Can't cancel" << std::endl;
#endif

      lMSB() = lf - EXTLONG_ONE;     // can't possibly cancel out
      sign() = sf;
    } else if (ls >= uf + EXTLONG_TWO) {
#ifdef CORE_TRACE
      std::cout << "Add/sub Rep:  Can't cancel" << std::endl;
#endif

      lMSB() = ls - EXTLONG_ONE;
      sign() = Op(ss);
    } else if (ffVal.isOK()) {// begin filter computation
#ifdef CORE_TRACE
      std::cout << "Add/sub Rep:  filter used" << std::endl;
#endif
#ifdef CORE_DEBUG_FILTER
      std::cout << "call filter in " << op() << "Rep" << std::endl;
#endif
      sign() = ffVal.sign();
      lMSB() = ffVal.lMSB();
      uMSB() = ffVal.uMSB();
    } else {                    // about the same size, might cancel out
#ifdef CORE_TRACE
      std::cout << "Add/sub Rep:  iteration start" << std::endl;
#endif

      extLong lowBound = computeBound();
      /* Zilin 06/11/2003
       * as BFMSS[2] might be a negative number, lowBound can be negative.
       * In this case, we just set it to 1 since we need at least one bit
       * to get the sign.  In the future, we may need to improve this.
       */
      if (lowBound <= EXTLONG_ZERO)
        lowBound = EXTLONG_ONE;

      if (!progressiveEvalFlag) {
        // convert the absolute error requirement "lowBound" to
        // a relative error requirement "ur", s.t.
        //    |x|*2^(-ur) <= 2^(-lowBound).
        // ==> r >= a + lg(x) >= a + (uMSB + 1);
        //          extLong  rf = lowBound + (uf + 1);
        //          extLong  rs = lowBound + (us + 1);
        //          first->approx(rf, CORE_INFTY);
        //          second->approx(rs, CORE_INFTY);
        // Chen: considering the uMSB is also an approximate bound.
        // we choose to use absolute precision up-front.
        Real newValue = Op(first->getAppValue(CORE_INFTY,
                                lowBound + EXTLONG_ONE),
                           second->getAppValue(CORE_INFTY,
                                   lowBound + EXTLONG_ONE));
        if (!newValue.isZeroIn()) { // Op(first, second) != 0
          lMSB() = newValue.lMSB();
          uMSB() = newValue.uMSB();   // chen: to get tighers value.
          sign() = newValue.sign();
        } else if (lowBound.isInfty()) {//check if rootbound is too big
          core_error("AddSubRep:root bound has exceeded the maximum size\n \
            but we still cannot decide zero.\n", __FILE__, __LINE__, false);
        } else {               // Op(first, second) == 0
          lMSB() = CORE_negInfty;
          sign() = 0;
        }
      } else {  // else do progressive evaluation
#ifdef CORE_TRACE
        std::cout << "Add/sub Rep:  progressive eval" << std::endl;
#endif
        // Oct 30, 2002: fixed a bug here!  Old versions used relative
        // precision bounds, but one should absolute precision for addition!
        // Moreover, this is much more efficient.

        // ua is the upper bound on the absolute precision in our iteration
        // Chee, Aug 8, 2004: it is important that ua be strictly
        //     larger than lowBound AND the defaultInitialProgressivePrec,
        //     so that we do at least one iteration of the for-loop. So:
        // i is the variable for iteration.
        extLong i = core_min(defInitialProgressivePrec, lowBound.asLong());
        extLong ua = lowBound.asLong() + EXTLONG_ONE;
        //   NOTE: ua is allowed to be CORE_INFTY
        
#ifdef CORE_DEBUG_BOUND
        std::cout << "DebugBound:" << "ua = " << ua << std::endl;
#endif
        // We initially set the lMSB and sign as if the value is zero:
        lMSB() = CORE_negInfty;
        sign() = 0;

        EscapePrecFlag = 0;     // reset the Escape Flag

        // Now we try to determine the real lMSB and sign,
        // in case it is not really zero:
#ifdef CORE_TRACE
        std::cout << "Upper bound (ua) for iteration is " << ua << std::endl;
        std::cout << "Starting iteration at i = " << i << std::endl;
#endif

        for ( ; i<ua; i*=EXTLONG_TWO) {
          // relative bits = i
          //
          // PROBLEM WITH NEXT LINE: you must ensure that
          // first and second are represented by BigFloats...
          //
          Real newValue = Op(first->getAppValue(CORE_INFTY, i),
                             second->getAppValue(CORE_INFTY, i));

#ifdef CORE_TRACE
          if (newValue.getRep().ID() == REAL_BIGFLOAT) 
          std::cout << "BigFloat! newValue->rep->ID() = "
                  << newValue.getRep().ID() << std::endl;
          else 
          std::cout << "ERROR, Not BigFloat! newValue->rep->ID() ="
                  << newValue.getRep().ID() << std::endl;
          std::cout << "newValue = Op(first,second) = "
                  << newValue << std::endl;
          std::cout << "first:appVal, appComputed, knownPrec, sign ="
                  << first->appValue() << ","
                  << first->appComputed() << ","
                  << first->knownPrecision() << ","
                  << first->sign() << std::endl;
          std::cout << "second:appVal, appComputed, knownPrec, sign ="
                  << second->appValue() << ","
                  << second->appComputed() << ","
                  << second->knownPrecision() << ","
                  << second->sign() << std::endl;
#endif
          if (!newValue.isZeroIn()) {   // Op(first, second) != 0
            lMSB() = newValue.lMSB();
            uMSB() = newValue.uMSB();
            sign() = newValue.sign();
#ifdef CORE_DEBUG_BOUND
            std::cout << "DebugBound(Exit Loop): " << "i=" << i << std::endl;
#endif
#ifdef CORE_TRACE
            std::cout << "Zero is not in, lMSB() = " << lMSB()
                    << ", uMSB() = " << uMSB()
                    << ", sign() = " << sign() << std::endl;
            std::cout << "newValue = " << newValue << std::endl;
#endif

            break; // assert -- this must happen in the loop if nonzero!
          }
          //8/9/01, Chee: implement escape precision here:
          if (i> EscapePrec) {
            EscapePrecFlag = -i.asLong();//negative means EscapePrec is used
            core_error("Escape precision triggered at",
                         __FILE__, __LINE__, false);
            if (EscapePrecWarning)
              std::cout<< "Escape Precision triggered at "
                      << EscapePrec << " bits" << std::endl;
#ifdef CORE_DEBUG
            std::cout << "EscapePrecFlags=" << EscapePrecFlag << std::endl;
            std::cout << "ua =" << ua  << ",lowBound=" << lowBound << std::endl;
#endif
            break;
          }// if
        }// for (long i=1...)

#ifdef CORE_DEBUG_BOUND
        rootBoundHitCounter++;
#endif

        if (sign() == 0 && ua .isInfty()) {
          core_error("AddSubRep: root bound has exceeded the maximum size\n \
            but we still cannot decide zero.\n", __FILE__, __LINE__, true);
        } // if (sign == 0 && ua .isInfty())
      }// else do progressive
    }
  }
  flagsComputed() = true;
}// AddSubRep::computeExactFlags

template <class Operator>
void AddSubRep<Operator>::computeApproxValue(const extLong& relPrec,
    const extLong& absPrec) {
  // Nov 13, 2002: added the analog of "reduceTo(first)" and "reduceTo(second)"
  //  that is found in computeExactFlags.  This is more efficient, but
  //  it also removes a NaN warning in subsequent logic!
  //  E.g., if first=0, then first->uMSB and first->lMSB are -infty, and
  //  subtracting them creates NaN.  Chee and Zilin.
  if (first->sign() == 0) {
    appValue() = Op(second->getAppValue(relPrec, absPrec));
    return;
  }
  if (second->sign() == 0) {
    appValue() = first->getAppValue(relPrec, absPrec);
    return;
  }
  if (lMSB() < EXTLONG_BIG && lMSB() > EXTLONG_SMALL) {
    extLong rf = first->uMSB()-lMSB()+relPrec+EXTLONG_FOUR;  // 2 better
    if (rf < EXTLONG_ZERO)
      rf = EXTLONG_ZERO;  // from Koji's thesis P63: Proposition 26
    extLong rs = second->uMSB()-lMSB()+relPrec+EXTLONG_FOUR; // 2 better
    if (rs < EXTLONG_ZERO)
      rs = EXTLONG_ZERO;  // from Koji's thesis P63: Proposition 26
    extLong  a  = absPrec + EXTLONG_THREE;                      // 1 better
    appValue() = Op(first->getAppValue(rf, a), second->getAppValue(rs, a));
  } else {
    std::cerr << "lMSB = " << lMSB() << std::endl; // should be in core_error
    core_error("CORE WARNING: a huge lMSB in AddSubRep",
                __FILE__, __LINE__, false);
  }
}

typedef AddSubRep<Add> AddRep;

typedef AddSubRep<Sub> SubRep;

class MultRep : public BinOpRep {
public:


  MultRep(ExprRep* f, ExprRep* s) : BinOpRep(f, s) {
    ffVal = first->ffVal * second->ffVal;
  }
  ~MultRep() {}

  CORE_MEMORY(MultRep)
protected:
  void computeExactFlags();
  void computeApproxValue(const extLong&, const extLong&);
  const std::string op() const {
    return "*";
  }
};

class DivRep : public BinOpRep {
public:


  DivRep(ExprRep* f, ExprRep* s) : BinOpRep(f, s) {
    ffVal = first->ffVal / second->ffVal;
  }
  ~DivRep() {}

  CORE_MEMORY(DivRep)
protected:
  void computeExactFlags();
  void computeApproxValue(const extLong&, const extLong&);
  const std::string op() const {
    return "/";
  }
};

// inline functions
inline int ExprRep::getExactSign() {
  if (!nodeInfo)
    initNodeInfo();

  if (!flagsComputed()) {
    degreeBound();
#ifdef CORE_DEBUG
    dagSize();
    fullClearFlag();
#endif
    computeExactFlags();
  }
  return sign();
}

// Chee, 7/17/02: degreeBound() function is now
// taken out of "computeExactFlags()
inline int ExprRep::getSign() {
  if (ffVal.isOK())
    return ffVal.sign();
  else
    return getExactSign();
}

// you need force to approximate before call these functions!!
inline BigInt ExprRep::BigIntValue() {
  return getAppValue().BigIntValue();
}

inline BigRat ExprRep::BigRatValue() {
  return getAppValue().BigRatValue();
}

inline BigFloat ExprRep::BigFloatValue() {
  return getAppValue().BigFloatValue();
}

CORE_END_NAMESPACE
#endif // _CORE_EXPRREP_H_