/****************************************************************** * Core Library Version 1.6, June 2003 * Copyright (c) 1995-2002 Exact Computation Project * * File: ExprRep.h * * Written by * Koji Ouchi * Chee Yap * Igor Pechtchanski * Vijay Karamcheti * Chen Li * Zilin Du * Sylvain Pion * Vikram Sharma * WWW URL: http://cs.nyu.edu/exact/ * Email: exact@cs.nyu.edu * * $Id$ *****************************************************************/ #ifndef CORE_EXPRREP_H #define CORE_EXPRREP_H #include #include #include #include #include #include CORE_BEGIN_NAMESPACE #ifdef 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 /// \struct NodeInfo /// \brief store information of a node struct NodeInfo { Real appValue; ///< current approximate value bool appComputed; ///< true if the approx value been computed bool flagsComputed; ///< true if rootBound parameters have been computed extLong knownPrecision; ///< Precision achieved by current approx value #ifdef DEBUG extLong relPrecision; extLong absPrecision; unsigned long numNodes; #endif /// d_e bounds the degree of the minimal polynomial of a DAG expression /** Basically, d_e is equal to 2^k where k is the number of square-root nodes * in the DAG. If there are other kinds of non-linear nodes, this is * generalized accordingly. */ extLong d_e; bool visited; ///< flag in counting # of sqrts int sign; ///< sign of the value being represented. extLong uMSB; ///< upper bound of the position of Most Significant Bit extLong lMSB; ///< lower bound of the position of Most Significant Bit // 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; ///< measure is really lg(Measure) // For our new bound. /// 2^{high(E)} is an UPPER bound for the moduli of ALL conjugates of E. /** In our papers, high is equal to log_2(\overline{\mu(E)}). */ extLong high; /// 2^{-low(E)} is an LOWER bound for the moduli of ALL NON_ZERO conjugate of E. /** BE CAREFUL! NOTE THAT UNLIKE "high", the sign of low is negated here! In our papers, low is equal to -log_2(\underline{\nu(E)}) */ extLong low; /// \brief upper bound of the leading coefficient of minimal defining /// polynomial of $E$. extLong lc; /// \brief upper bound of the last non-zero coefficient of minimal defining /// polynomial of $E$. extLong tc; // For the 2-ary BFMSS bound. extLong v2p, v2m; // For the 5-ary BFMSS bound. extLong v5p, v5m; /// 2^u25 is an upper bound for the moduli of all the conjugates of U(E) /** where E = 2^v2*5^v5*U(E)/L(E), U(E) and L(E) are division-free. */ extLong u25; /// 2^l25 is an upper bound for the moduli of all the conjugates of L(E) /** where E = 2^v2*5^v5*U(E)/L(E), U(E) and L(E) are division-free. */ extLong l25; int ratFlag; ///< rational flag BigRat* ratValue; ///< rational value /// default constructor NodeInfo(); }; // forward reference // class Expr; /// \class ExprRep /// \brief The sharable, internal representation of expressions class ExprRep { public: /// \name Constructor and Destructor //@{ /// default constructor ExprRep(); /// virtual destructor for this base class virtual ~ExprRep() { if (nodeInfo != NULL) // This check is only for optimization. delete nodeInfo; } //@} /// \name Reference Counting //@{ /// increase reference counter void incRefCount() { ++refCount; } /// decrease reference counter void decRefCount() { if ((--refCount) == 0) delete this; } /// check whether reference counter == 1 int isUnique() const { return refCount == 1; } //@} /// \name Helper Functions //@{ /// Get the approximate value const Real & getAppValue(const extLong& relPrec = defRelPrec, const extLong& absPrec = defAbsPrec); /// Get the sign. 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 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; } /// Get BigFloat BigInt BigIntValue(); BigRat BigRatValue(); BigFloat BigFloatValue(); /// represent as a string in decimal value // toString() Joaquin Grech 31/5/2003 std::string toString(long prec=defOutputDigits, bool sci=false) const { return nodeInfo->appValue.toString(prec,sci); } //@} /// \name Debug functions //@{ /// dump the contents in this DAG node const std::string dump(int = OPERATOR_VALUE) const; /// print debug information in list mode virtual void debugList(int level, int depthLimit) const = 0; /// print debug information in tree mode virtual void debugTree(int level, int indent, int depthLimit) const = 0; //@} /// \name I/O Stream //@{ friend std::ostream& operator<<(std::ostream&, ExprRep&); //@} CORE_MEMORY(ExprRep) private: unsigned refCount; // reference count public: enum {OPERATOR_ONLY, VALUE_ONLY, OPERATOR_VALUE, FULL_DUMP}; NodeInfo* nodeInfo; ///< node information filteredFp ffVal; ///< filtered value /// \name Approximation Functions //@{ /// initialize nodeInfo virtual void initNodeInfo() = 0; /// compute the sign, uMSB, lMSB, etc. virtual void computeExactFlags() = 0; /// compute the minimal root bound extLong computeBound(); /// driver function to approximate void approx(const extLong& relPrec, const extLong& absPrec); /// compute an approximate value satifying the specified precisions virtual void computeApproxValue(const extLong&, const extLong&) = 0; /// Test whether the current approx. value satisfies [relPrec, absPrec] bool withinKnownPrecision(const extLong&, const extLong&); //@} /// \name Misc Functions //@{ /// reduce current node 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 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 ~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; //@} CORE_MEMORY(ConstRep) 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(); } /// clear visited flag void clearFlag() { visited() = false; } #ifdef 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 ConstPolyRep : public ConstRep { public: /// \name Constructors and Destructor //@{ /// default constructor ConstPolyRep() { } /// constructor for Polynomial ConstPolyRep(const Polynomial& 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) { std::cerr << "error! root index out of bound." << std::endl; abort(); } // test if the root isolated in I is 0: if ((I.first == 0)&&(I.second == 0)) ffVal = 0; else ffVal = computeFilteredValue(); } /// constructor for Polynomial ConstPolyRep(const Polynomial& p, const BFInterval& II) : ss(p), I(II) { if (ss.numberOfRoots(I.first, I.second) != 1) { std::cerr << "error! non-isolating interval." << std::endl; abort(); } ffVal = computeFilteredValue(); } /// destructor ~ConstPolyRep() {} //@} CORE_MEMORY(ConstPolyRep) private: Sturm ss; ///< internal Sturm sequences BFInterval I; ///< current interval contains the real value 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 here). // Hence, the precision we pass to refine should be 54-lMSB(k). // refine with newton (new method) I = ss.newtonRefine(I, 54-(ss.seq[0].CauchyLowerBound()).lMSB().asLong()); //return I.first.doubleValue(); // NOTE: This is not quite right! // It should be "centralize" which should set // the error bit correctly. // E.g., 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(); 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 = ((BigInt(x.err()) << 53) / (x.m() + x.err())).longValue(); return filteredFp(val, max, ind); } 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!! // 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->incRefCount(); } /// destructor virtual ~UnaryOpRep() { child->decRefCount(); } //@} /// \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; //@} CORE_MEMORY(UnaryOpRep) protected: ExprRep* child; ///< pointer to its child node /// initialize nodeInfo virtual void initNodeInfo(); /// clear visited flag void clearFlag(); #ifdef 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->incRefCount(); second->incRefCount(); } /// destructor virtual ~BinOpRep() { first->decRefCount(); second->decRefCount(); } //@} /// \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; //@} CORE_MEMORY(BinOpRep) 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 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 const T& operator()(const T& t) const { return t; } /// binary operator template 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 T operator()(const T& t) const { return -t; } /// binary operator template 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 AddSubRep : public BinOpRep { public: /// \name Constructors and Destructor //@{ /// constructor AddSubRep(ExprRep* f, ExprRep* s) : BinOpRep(f, s) { ffVal = Op(first->ffVal, second->ffVal); } /// destructor ~AddSubRep() {} //@} CORE_MEMORY(AddSubRep) 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 Operator::name; } private: static Operator Op; }; template Operator AddSubRep::Op; /// \typedef AddRep /// \brief AddRep for easy of use typedef AddSubRep AddRep; /// \typedef SubRep /// \brief SuRep for easy of use typedef AddSubRep SubRep; /// \class MultRep /// \brief multiplication operator node class MultRep : public BinOpRep { public: /// \name Constructors and Destructor //@{ /// constructor MultRep(ExprRep* f, ExprRep* s) : BinOpRep(f, s) { ffVal = first->ffVal * second->ffVal; } /// destructor ~MultRep() {} //@} CORE_MEMORY(MultRep) 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 "*"; } }; /// \class DivRep /// \brief division operator node class DivRep : public BinOpRep { public: /// \name Constructors and Destructor //@{ /// constructor DivRep(ExprRep* f, ExprRep* s) : BinOpRep(f, s) { ffVal = first->ffVal / second->ffVal; } /// destructor ~DivRep() {} //@} CORE_MEMORY(DivRep) 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 "/"; } }; // inline functions inline int ExprRep::getExactSign() { if (!nodeInfo) initNodeInfo(); if (!flagsComputed()) { degreeBound(); #ifdef 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