upgrade to v1.6

Packages/Core/include/CORE/geom2d/circle2d.h

@@ -0,0 +1,132 @@
+/*****************************************************************
+ * File: circle2d.h
+ * Synopsis:
+ *      Basic 2-dimensional geometry
+ * Author: Shubin Zhao (shubinz@cs.nyu.edu), 2001.
+ *
+ *****************************************************************
+ * CORE Library Version 1.4 (July 2001)
+ *       Chee Yap <yap@cs.nyu.edu>
+ *       Chen Li <chenli@cs.nyu.edu>
+ *       Zilin Du <zilin@cs.nyu.edu>
+ *
+ * Copyright (c) 1995, 1996, 1998, 1999, 2000, 2001 Exact Computation Project
+ *
+ * WWW URL: http://cs.nyu.edu/exact/
+ * Email: exact@cs.nyu.edu
+ *
+ * $Id$
+ *****************************************************************/
+
+#ifndef _CIRCLE2D_H
+#define _CIRCLE2D_H
+
+#include "point2d.h"
+#include "line2d.h"
+
+class Circle2d : public GeomObj {
+
+   /* An instance C of the data type circle is an oriented circle 
+      in the plane passing through three points p1, p2, p3. The 
+      orientation of C is equal to the orientation of the three defining
+      points. i.e. orientation(p1, p2, p3).
+      If \Labs{\{p1, p2, p3\}} = 1, C is the empty circle with center p1.
+      If p1, p2 and p3 are collinear, C is a straight line passing through
+      p1, p2 and p3 in this order and the center of C is undefined.    
+    */
+
+private:
+ 
+   Point2d p1;  // the 3 points defining the circle
+   Point2d p2;
+   Point2d p3;
+
+   int orient;  //orientation(p1, p2, p3)
+
+   Point2d* cp;  //pointer to center
+   double * rp;  //pointer to radius
+
+public:
+
+   Circle2d( const Point2d& p1, const Point2d& p2,  const Point2d& p3);
+   //initialized to the oriented circle through points p1, p2, p3 
+
+   Circle2d(const Point2d& a, const Point2d& b0);
+   //initialized to the counter-clockwise oriented circle with center a 
+   //passing through b0 
+
+   Circle2d(const Point2d& p);
+   //initialized to the trivial circle with center p
+
+   Circle2d();
+   //initialized to the trivial circle with center (0,0)
+
+   Circle2d(const Point2d& c, double r);
+   //initialized to the circle with center c and radius r with positive
+   //(i.e. counter-clockwise) orientation
+
+   Circle2d(const Circle2d& c);
+   //copy constructor 
+
+   virtual ~Circle2d(); 
+
+   Circle2d& operator=(const Circle2d& C);
+
+   //operations
+
+   Point2d center();
+   //return the center of the circle
+
+   double radius();
+   //returns the radius. 
+   //precond: the orientation of the circle is not 0
+
+   Point2d point1() const { return p1; }
+   Point2d point2() const { return p2; }
+   Point2d point3() const { return p3; }
+
+//   Point2d point_on_circle(float alpha);
+   //returns a point p on the circle with angle of alpha
+
+   bool is_degerate() const { return orient == 0; }
+   //returns true if the defining points are collinear
+
+   bool is_trivial() const {return p1 == p2; }
+   //returns true if radius is zero
+
+   int orientation() const { return orient; }
+
+   int side_of(const Point2d& p) const; 
+   // returns -1, +1 or 0 if p lies right of, left of or on the circle
+   // respectively
+
+   bool inside(const Point2d& p);
+   //returns true if p lies inside of the circle
+
+   bool outside(const Point2d& p); 
+   
+   bool contains(const Point2d& p) const ;
+   //returns true if p lies on the circle, false otherwise
+
+   double distance(const Point2d& p);
+   //returns the distance between p and the circle: distance to center - radius
+
+   double distance(const Line2d& l);
+   //returns the distance between l and the circle
+   //distance from center to l minus radius
+
+   double distance(Circle2d& D);
+   //returns the distance between this circle and circle D
+   //distance between two centers minus two radius
+
+   bool operator==(const Circle2d& D) const ;
+   bool operator!=(const Circle2d& D) const 
+	{ return !operator==(D); }
+
+   friend std::ostream& operator<<(std::ostream& out, Circle2d& c);
+   friend std::istream& operator>>(std::istream& in, Circle2d c); //?? Circle2d &
+
+}; // class Circle2d
+
+
+#endif

Packages/Core/include/CORE/geom2d/line2d.h

@@ -0,0 +1,152 @@
+/*****************************************************************
+ * File: line2d.h
+ * Synopsis:
+ *      Basic 2-dimensional geometry
+ * Author: Shubin Zhao (shubinz@cs.nyu.edu), 2001.
+ *
+ *****************************************************************
+ * CORE Library Version 1.4 (July 2001)
+ *       Chee Yap <yap@cs.nyu.edu>
+ *       Chen Li <chenli@cs.nyu.edu>
+ *       Zilin Du <zilin@cs.nyu.edu>
+ *
+ * Copyright (c) 1995, 1996, 1998, 1999, 2000, 2001 Exact Computation Project
+ *
+ * WWW URL: http://cs.nyu.edu/exact/
+ * Email: exact@cs.nyu.edu
+ *
+ * $Id$
+ *****************************************************************/
+
+#ifndef _LINE2D_H_
+#define _LINE2D_H_
+
+#include "point2d.h"
+
+class Line2d : public GeomObj {
+
+  /* An instance l of the data type $line$ is a directed straight line
+     in the two dimensional plane. The angle between a right oriented 
+     horizontal line and $l$ is called the direction of $l$.
+   */
+  /* member Vector is not used in this class, it's intended for use in
+       the operator+,- etc 
+     need to do: assure p0 != p1
+   */
+private:
+
+  Point2d p0;
+  Point2d p1;
+  Vector V;
+
+public:
+
+  /*************************************************************
+   *  constructors
+   *************************************************************/
+
+  Line2d(const Point2d & p, const Vector &v);
+  // line initialized to pass through points p and p+v
+
+  Line2d(const Point2d &p, const Point2d &q);
+  //line is initialized to pass through points p and q directed from p to q
+
+//  Line2d(const point& p, double alpha);
+  //line passes through point p with direction alpha
+
+  Line2d(const Line2d &);
+
+  Line2d();
+  //line passes through the origin with direction 0.
+
+  virtual ~Line2d() {}
+  /*************************************************************
+   *   member functions
+   *************************************************************/
+
+  Vector direction() const { return p1-p0; }
+   // returns the direction as a vector 
+
+  Point2d startPt() const { return p0; }  
+  Point2d stopPt() const  { return p1; }
+
+  double distance(Point2d q) const;
+   // returns the Euclidean distance between this line and point q
+
+  Point2d projection(const Point2d& p) const;
+   // returns the projection of p on this line
+
+  int orientation( const Point2d& p ) const;
+   // orientation of p0, p1 and p
+
+   // the sine/cosine of the angle made with positive x-direction
+  double sine() const { return (p1.Y() - p0.Y()) / p0.distance(p1); }
+  double cosine() const { return (p1.X() - p0.X()) / p0.distance(p1); }
+
+  Line2d rotate90( const Point2d& q)
+  { return Line2d(startPt().rotate90(q), stopPt().rotate90(q)); }
+
+  double y_abs() const;
+  // returns the y-abscissa of the line 
+
+  double slope() const ;    
+  //precond: is not vertical
+
+  /*************************************************************
+   *   predicates
+   *************************************************************/
+
+  bool isVertical() const { return p0.X() == p1.X(); }
+  bool isHorizontal() const { return p0.Y() == p1.Y(); }
+  bool is_trivial() const {return p0 == p1; }   //meaning for a line?
+
+  bool contains( const Point2d& p) const { 
+          return orientation2d(p0, p1, p) == 0; }
+  bool isCoincident( const Line2d& g) const { 
+          return contains(g.p0) && contains(g.p1); }  
+  bool isParallel(const Line2d& l) const {
+    return det(V, l.direction()) == 0; }
+
+  bool operator==( const Line2d& g ) const { return isCoincident(g); }
+  bool operator!=( const Line2d& g ) const { return !operator==(g); }
+
+  /*************************************************************
+   *   intersection
+   *************************************************************/
+
+  int intersects(const Line2d& t) const;
+   // decides whether *this and t intersects
+   // return dim of intersection. 
+   // return -1 if no intersection
+
+  GeomObj* intersection(const Line2d& g) const;
+   //if this line and g intersect in a single point, this point is 
+   // assigned to p and the result is true, otherwise the result is false
+
+  /*************************************************************
+   *   angles and others
+   *************************************************************/
+  friend int orientation2d( const Line2d& l, const Point2d& p);
+  // computes the orientation (a, b, p), where a!=b and a and b appear 
+  // in this order on line l
+
+  friend int cmp_slopes(const Line2d& l1, const Line2d& l2) 
+  //l1.slope > l2.slope: +1; equal: 0; otherwise: -1
+  {
+     if (l1.slope() == l2.slope())
+         return 0;
+     else
+         return (l1.slope() > l2.slope()) ? +1 : -1;
+  } 
+  
+  /*************************************************************
+   *   I/O
+   *************************************************************/
+
+  friend std::istream& operator>>(std::istream& in, Line2d& l);
+  friend std::ostream &operator<<(std::ostream & out, const Line2d & l);
+}; // class Line2d
+
+extern Line2d p_bisector(Point2d& p, Point2d& q);
+
+#endif

Packages/Core/include/CORE/geom2d/point2d.h

@@ -0,0 +1,89 @@
+/*****************************************************************
+ * File: point2d.h
+ * Synopsis:
+ *      Basic 2-dimensional geometry
+ * Author: Shubin Zhao (shubinz@cs.nyu.edu), 2001.
+ *
+ *****************************************************************
+ * CORE Library Version 1.4 (July 2001)
+ *       Chee Yap <yap@cs.nyu.edu>
+ *       Chen Li <chenli@cs.nyu.edu>
+ *       Zilin Du <zilin@cs.nyu.edu>
+ *
+ * Copyright (c) 1995, 1996, 1998, 1999, 2000, 2001 Exact Computation Project
+ *
+ * WWW URL: http://cs.nyu.edu/exact/
+ * Email: exact@cs.nyu.edu
+ *
+ * $Id$
+ *****************************************************************/
+
+
+#ifndef _POINT2D_H
+#define _POINT2D_H
+
+#ifndef CORE_LEVEL
+#  define CORE_LEVEL 3
+#endif
+
+#include "CORE.h"
+#include "../linearAlgebra.h"
+#include "../geombase.h"
+
+class Point2d : public GeomObj {
+
+private:
+  double x, y;
+
+public:
+ 
+  //constructors
+  Point2d();  //initialized to origin(0,0)
+  Point2d(double, double);
+  Point2d(const Point2d &);
+  Point2d(Vector v);
+  //create a point initialized to the point $(v[0], v[1])$ 
+  //precondition: v.dim() = 2
+
+  //destrutor
+  virtual ~Point2d() {}
+
+  Point2d& operator=(const Point2d&);
+  
+  double X() const { return x; }
+  double Y() const { return y; } 
+
+  Vector toVector() const { return Vector(X(), Y()); } 
+
+  int dim() const { return 2; }
+
+  double distance(const Point2d) const;
+  // returns the Euclidean distance between p and this
+
+  double distance() const { return distance(Point2d(0, 0)); }
+  // returns distance between this and origin
+
+  Vector operator-(const Point2d &) const;
+  Point2d operator+(const Vector &) const;
+
+  Point2d rotate90( const Point2d& q);
+  // returns the point rotated about q by angle of 90 degrees
+
+  bool operator==(const Point2d&) const;
+  bool operator!=(const Point2d& p) const {return !operator==(p); }
+  
+  friend std::ostream& operator<< (std::ostream&, const Point2d);
+  // write point p to output stream
+
+  friend std::istream& operator>>(std::istream&, Point2d&);
+  // reads the x and y coordinates of point p from the input stream
+
+};
+
+
+Point2d center(Point2d& a, Point2d& b); 
+int orientation2d(const Point2d& a, const Point2d& b, const Point2d& c);
+double area(const Point2d& a, const Point2d& b, const Point2d& c);
+bool collinear(const Point2d& a, const Point2d& b, const Point2d& c);
+
+#endif

Packages/Core/include/CORE/geom2d/segment2d.h

@@ -0,0 +1,177 @@
+/*****************************************************************
+ * File: segment2d.h
+ * Synopsis:
+ *      Basic 2-dimensional geometry
+ * Author: Shubin Zhao (shubinz@cs.nyu.edu), 2001.
+ *
+ *****************************************************************
+ * CORE Library Version 1.4 (July 2001)
+ *       Chee Yap <yap@cs.nyu.edu>
+ *       Chen Li <chenli@cs.nyu.edu>
+ *       Zilin Du <zilin@cs.nyu.edu>
+ *
+ * Copyright (c) 1995, 1996, 1998, 1999, 2000, 2001 Exact Computation Project
+ *
+ * WWW URL: http://cs.nyu.edu/exact/
+ * Email: exact@cs.nyu.edu
+ *
+ * $Id$
+ *****************************************************************/
+
+
+#ifndef _SEGMENT2D_H_
+#define _SEGMENT2D_H_
+
+#include "point2d.h"
+#include "line2d.h"
+
+/************************************************************
+ *  Class Segment2d:
+ *
+ *  An instance s of Segment2d is a finite or infinite line segment
+ *  in the two dimensional plane, defined by a start point
+ *  s.startPt() and a stop point s.stopPt().  It can be regarded
+ *  as an open or a closed segment (default: open), and directed
+ *  or not (default: directed).
+ *
+ *  We do not necessarily assume that startPt() != stopPt().
+ ************************************************************/
+
+class Segment2d : public GeomObj {
+
+private:
+
+  Point2d p0;
+  Point2d p1;
+  bool directed;  // segments can be directed or not (default is directed)
+  bool open;     // segments can be open or closed (default is open)
+
+public:
+
+  /************************************************************
+   *   constructors
+   ************************************************************/
+
+  Segment2d(const Segment2d &);
+
+  Segment2d(const Point2d &p, const Point2d &q);
+  	//finite segment with endpoints p and q
+
+  Segment2d(const Point2d & p, const Vector & v);
+	//ray segment
+
+  Segment2d();
+  	//unit segment from (0,0) to (1,0) 
+
+  virtual ~Segment2d() {}
+  /*************************************************************
+   *   member functions
+   *************************************************************/
+
+  Point2d startPt() const { return p0; }  
+  Point2d stopPt() const { return p1; }
+
+  void reverse() { Point2d pTmp = p0; p0=p1; p1=pTmp; }
+   //  reverses the direction of the segment
+
+  Line2d toLine() const { return Line2d(p0,p1); }
+ 
+  double length() const { return p0.distance(p1); }
+    //length of segment
+
+  double distance( const Point2d& p ) const;
+   // returns the Euclidean distance between this segment and point q
+
+  Point2d nearPt( const Point2d& p ) const;
+   // returns the point on segment closest to q;
+
+  void setDirected( bool _directed ) { directed = _directed; }
+  void setOpen( bool _open ) { directed = open; }
+
+   // orientation of p0, p1 and p
+  int orientation( const Point2d& p ) const {
+    return toLine().orientation(p); }
+
+  /*************************************************************
+   *   predicates
+   *************************************************************/
+  bool isOpen() const {return open; }
+  bool isDirected() const {return directed; }
+  bool isTrivial() const {return p0 == p1; }  
+  bool isVertical() const { return p0.X() == p1.X(); }
+  bool isHorizontal() const { return p0.Y() == p1.Y(); }
+  bool isCollinear( Point2d& p ) const { return toLine().contains(p); }
+  bool isCoincident( const Segment2d& s) const;
+  bool isParallel( const Segment2d& s ) {
+    return toLine().isParallel( s.toLine() );  }
+
+  bool contains( const Point2d& p ) const;
+  bool contains( const Segment2d& s ) const { 
+     return contains(s.startPt()) && contains(s.stopPt()); }  
+
+  bool operator==(const Segment2d& s) const { return isCoincident(s); }
+
+  bool operator!=(const Segment2d& s) const { return !isCoincident(s); }
+
+  /*************************************************************
+   *   intersection
+   *************************************************************/
+
+  int intersects( const Line2d& l ) const;
+   //decides whether *this and t intersect in one point
+   // return dim of intersetion
+  
+  int intersects( const Segment2d& s ) const;
+   //decides whether *this and t intersect in one point
+   // return dim of intersetion
+
+  GeomObj* intersection( const Line2d& l ) const;
+   // return intersection point if this segment and l intersect at a single point
+   // the intersection point is returned 
+  
+  GeomObj* intersection( const Segment2d& s ) const;
+   // return intersection point if this segment and s intersect at a single point
+   // the intersection point is returned 
+ 
+  /*************************************************************
+   *   angles
+   *************************************************************/
+
+   // the sine/cosine of the angle made with positive x-direction
+  double sine() const { return (p1.Y() - p0.Y()) / length() ; }
+  double cosine() const { return (p1.X() - p0.X()) / length() ; }
+
+  Line2d rotate90(const Point2d& q)
+  { return Line2d(startPt().rotate90(q), stopPt().rotate90(q)); }
+
+   // computes the orientation (a, b, p), where a!=b and a and b appear 
+   // in this order on segment l
+  friend int orientation2d( const Segment2d& s, const Point2d& p) {
+    return orientation2d( s.toLine(), p ); 
+  }
+
+  friend int cmp_slopes( const Segment2d& s1, const Segment2d& s2) 
+   //l1.slope > l2.slope: +1; equal: 0; otherwise: -1
+  {
+     Line2d l1 = s1.toLine();
+     Line2d l2 = s2.toLine();
+     if (l1.slope() == l2.slope())
+         return 0;
+     else
+         return (l1.slope() > l2.slope()) ? +1 : -1;
+  } 
+  
+  /*************************************************************
+   *   I/O
+   *************************************************************/
+
+
+  friend std::istream& operator>>(std::istream& in, Segment2d& l);
+  friend std::ostream &operator<<(std::ostream & out, const Segment2d & l);
+   // syntax: {[} p {===} q {]}
+
+}; //class Segment2d
+
+extern Line2d p_bisector(Point2d& p, Point2d& q);
+
+#endif

Packages/Core/include/CORE/poly/Poly.h

@@ -0,0 +1,327 @@
+/******************************************************************
+ * Core Library Version 1.6, June 2003
+ * Copyright (c) 1995-2003 Exact Computation Project
+
+ File: Poly.h
+
+ Description: simple polynomial class
+
+	REPRESENTATION:
+	--Each polynomial has a nominal "degree" (this
+		is an upper bound on the true degree, which
+		is determined by the first non-zero coefficient).
+	--coefficients are parametrized by some number type "NT".
+	--coefficients are stored in the "coeff" array of
+		length "degree + 1".  
+		CONVENTION: coeff[i] is the coefficient of X^i.  So, a
+			    coefficient list begins with the constant term.
+	--IMPORTANT CONVENTION:
+		the zero polynomial has degree -1
+		while nonzero constant polynomials have degree 0.
+
+	FUNCTIONALITY:
+	--Polynomial Ring Operations (+,-,*)
+	--Power
+	--Evaluation
+	--Differentiation
+	--Remainder, Quotient 
+	--Resultant, Discriminant (planned)
+	--Polynomial Composition (planned)
+	--file I/O (planned)
+	
+ Author: Chee Yap 
+ Date:   May 28, 2002
+ Core Library 1.4.1
+ $Id$
+ ************************************** */ 
+
+#ifndef CORE_POLY_H
+#define CORE_POLY_H
+
+#include "../BigFloat.h"
+#include <vector>
+
+CORE_BEGIN_NAMESPACE
+
+class Expr;
+// ==================================================
+// Typedefs 
+// ==================================================
+
+//typedef std::vector<Expr>	VecExpr;
+//typedef std::pair<Expr, Expr>	Interval;
+//typedef std::vector<Interval>	VecInterval;
+typedef std::pair<BigFloat, BigFloat>	BFInterval;
+			// NOTE: an error condition is indicated by
+			// the special interval (1, 0)
+typedef std::vector<BFInterval>	BFVecInterval;
+
+
+// ==================================================
+// Polynomial Class
+// ==================================================
+
+template <class NT>
+class Polynomial {
+public:
+  typedef std::vector<NT> VecNT;
+
+  int degree;	// This is the nominal degree (an upper bound 
+		// on the true degree)
+  NT * coeff;	// coeff is an array of size degree+1;
+		//	This remark holds even when degree = -1.
+		// Notes:
+		// (1) coeff[i] is the coefficient of x^i
+		// (2) The Zero Polynomial has degree -1
+		// (3) Nonzero Constant Polynomials has degree 0
+
+  // STATIC MEMBERS
+  // static NT ccc_; // THIS IS A TEMPORARY HACK
+  static  int COEFF_PER_LINE;		// pretty print parameters
+  static const char * INDENT_SPACE;		// pretty print parameters
+
+  static const Polynomial<NT> & polyZero();
+  static const Polynomial<NT> & polyUnity();
+  static Polynomial polyWilkinson;	      // a sample polynomial
+
+  // Constructors:
+  Polynomial(void);	// the Zero Polynomial
+  Polynomial(int n);	// construct the Unit Polynomial of nominal deg n>=0
+  Polynomial(int n, NT * coef);
+  Polynomial(const Polynomial &);
+  Polynomial(const VecNT &);
+  Polynomial(int n, const char* s[]);
+  ~Polynomial();
+
+  // Assignment:
+  Polynomial & operator=(const Polynomial&);
+
+  // Expand and Contract  
+  //  -- they are semi-inverses: i.e., Contract(expand(p))=p
+  int expand(int n);	// Change the current degree to n
+			// Helper function for polynomial arithmetic
+  int contract();		// get rid of leading zeros
+
+  // Polynomial arithmetic (these are all self-modifying):
+  Polynomial & operator+=(const Polynomial&);	// +=
+  Polynomial & operator-=(const Polynomial&);	// -=
+  Polynomial & operator*=(const Polynomial&);	// *=
+  Polynomial & operator-();			// unary minus
+  Polynomial & power (unsigned int n) ;		// power (*this is changed!)
+
+  Polynomial & mulScalar (const NT & c);	// return (*this) * (c)
+  Polynomial & mulXpower(int i); 	// if i >= 0, then this is equivalent
+					//	to multiplying by X^i
+					// if i < 0 to dividing by X^i
+  Polynomial pseudoRemainder (Polynomial& p);
+			// the pseudo quotient of (*this) mod p
+			//	is returned, but (*this) is transformed
+			//	into the pseudo remainder
+  Polynomial reduceStep (Polynomial& p );
+			// One step of pseudo remainder 
+			// What is returned is a special polynomial C+M
+			//	telling us the initial constant C and					//	the quotient M of C*(THIS) divided by p.
+			
+  // Polynomial gcd(Polynomial p); // This may not be defined for some NT domains
+
+  // Get functions
+  int getDegree() const;	// nominal degree
+  int getTrueDegree() const;	// true degree
+  const NT & getLeadCoeff() const;	// get TRUE leading coefficient
+  const NT & getTailCoeff() const;	// get last non-zero coefficient
+  NT** getCoeffs() ;		// get all coefficients
+  const NT& getCoeff(int i) const;	// Get single coefficient of X^i
+				// NULL pointer if invalid i
+  // Set functions
+  bool setCoeff(int i, const NT & cc);	// Make cc the coefficient of X^i
+					// Return FALSE if invalid i
+					// !! User's responsibility to 
+					// delete the old coefficient if
+					// necessary !!
+  // Helper Functions
+  void reverse();		// reverse the coefficients;  useful when
+				// input coefficients are in reversed
+
+  // Evaluation
+  Expr eval(const Expr&) const;		// evaluation
+  BigFloat eval(const BigFloat&) const;	// evaluation
+
+  // Bounds
+  BigFloat CauchyUpperBound() const;  // Cauchy Root Upper Bound
+  BigFloat CauchyLowerBound() const;  // Cauchy Root Lower Bound
+  BigFloat sepBound() const;	// separation bound (multiple roots allowed)
+  BigFloat height() const;	// height return type BigFloat
+  BigFloat length() const;	// length return type BigFloat
+
+  // Differentiation
+  Polynomial & differentiate() ;		// self-differentiation
+  Polynomial & differentiate(int n) ;		// multi self-differentiation
+
+  // Reductions of polynomials
+  Polynomial & squareFreePart(); // P/gcd(P,P')
+  Polynomial & primPart(); 	//  Primitive Part of *this (which is changed)
+
+  //////////////////////////////////////////////////////////////////
+  // Resultant and discriminant
+  // NT & resultant() ;				// resultant
+  // NT & discriminant() ;				// discriminant
+
+  // Composition of Polynomials:
+  // NOT yet implemented
+
+  //////////////////////////////////////////////////////////////////
+  // Polynomial Dump
+  void dump() const;			// plain dump
+  void dump(std::string m1) const;	// dump with message string
+  void mapleDump() const;		// dump of maple code for Polynomial
+
+}; //Polynomial Class
+
+// template < class NT >
+//     NT Polynomial<NT>::ccc_;
+
+// ==================================================
+// Static Constants
+//	Does this belong here?
+// ==================================================
+
+template < class NT >
+  CORE_INLINE
+  const Polynomial<NT> & Polynomial<NT>::polyZero(){
+	static Polynomial<NT> zeroP;
+	return zeroP;
+  }
+
+template < class NT >
+  CORE_INLINE
+  const Polynomial<NT> & Polynomial<NT>::polyUnity() {
+	static NT c[] = {1};
+	static Polynomial<NT> unityP(0, c);
+	return unityP;
+  }
+
+// ==================================================
+// Useful functions for Polynomial class
+// ==================================================
+
+  // polynomial arithmetic:
+template < class NT >
+  Polynomial<NT> operator+(const Polynomial<NT>&, const Polynomial<NT>&);// +
+template < class NT >
+  Polynomial<NT> operator-(const Polynomial<NT>&, const Polynomial<NT>&);// -
+template < class NT >
+  Polynomial<NT> operator*(const Polynomial<NT>&, const Polynomial<NT>&);// *
+template < class NT >
+  Polynomial<NT> power(const Polynomial<NT>&, int n);		// power
+template < class NT >
+  Polynomial<NT> differentiate(const Polynomial<NT>&);		// differentiate
+template < class NT >
+  Polynomial<NT> differentiate(const Polynomial<NT>&, int n);	// multi-differ.
+
+  // comparisons
+template < class NT >
+  bool operator==(const Polynomial<NT>&, const Polynomial<NT>&); // ==
+template < class NT >
+  bool operator!=(const Polynomial<NT>&, const Polynomial<NT>&); // !=
+template < class NT >
+  bool zeroP(const Polynomial <NT>&);			// =Zero Poly?
+template < class NT >
+  bool unitP(const Polynomial <NT>&);			// =Unit Poly?
+
+  // stream i/o
+template < class NT >
+  std::ostream& operator<<(std::ostream&, const Polynomial<NT>&);
+template < class NT >
+  std::istream& operator>>(std::istream&, Polynomial<NT>&);
+
+// ==================================================
+// Inline Functions
+// ==================================================
+
+  // friend polynomial arithmetic:
+template < class NT >
+  CORE_INLINE
+  Polynomial<NT> operator+(const Polynomial<NT>& p1,
+				const Polynomial<NT>& p2){	// +
+    return Polynomial<NT>(p1) += p2;
+  }
+template < class NT >
+  CORE_INLINE
+  Polynomial<NT> operator-(const Polynomial<NT>& p1,
+				const Polynomial<NT>& p2){	// -
+    return Polynomial<NT>(p1) -= p2;
+  }
+template < class NT >
+  CORE_INLINE
+  Polynomial<NT> operator*(const Polynomial<NT>& p1,
+				const Polynomial<NT>& p2){	// *
+    return Polynomial<NT> (p1) *= p2;
+  }
+template < class NT >
+  CORE_INLINE
+  Polynomial<NT> power(const Polynomial<NT>& p, int n){				// power
+    return Polynomial<NT>(p).power(n);
+  }
+
+// equal to zero poly?
+template < class NT >
+  CORE_INLINE
+  bool zeroP(const Polynomial <NT>& p){			// =Zero Poly?
+	return (p.getTrueDegree()== -1);
+  }
+template < class NT >
+  CORE_INLINE
+  bool unitP(const Polynomial <NT>& p){			// =Unit Poly?
+	int d = p.getTrueDegree();
+	return ((d == 0) && p.coeff[0]==1 );
+  } 
+
+// get functions
+template < class NT >
+  CORE_INLINE
+  int Polynomial<NT>::getDegree() const {
+	return degree;
+  }
+// get TRUE leading coefficient
+template < class NT >
+  CORE_INLINE
+  const NT & Polynomial<NT>::getLeadCoeff() const {
+	return getCoeff(getTrueDegree());
+  }
+
+// get last non-zero coefficient
+template < class NT >
+  CORE_INLINE
+  const NT & Polynomial<NT>::getTailCoeff() const {
+     for (int i = 0; i<= getTrueDegree(); i++) 
+	if (coeff[i] != 0) return coeff[i];
+     // This ought to be an error (user should check this :
+     NT * zero = new NT(0);
+     return *zero;
+  }
+
+template < class NT >
+  CORE_INLINE
+  NT** Polynomial<NT>::getCoeffs() {
+	return &coeff;
+  }
+template < class NT >
+  CORE_INLINE
+  const NT& Polynomial<NT>::getCoeff(int i) const {
+	//if (i > degree) return NULL;
+	assert(i <= degree);
+	return coeff[i];
+  }
+// set functions
+template < class NT >
+  CORE_INLINE
+  bool Polynomial<NT>::setCoeff(int i, const NT & cc) {
+	if ((i<0) || (i > degree)) return false;
+	coeff[i] = cc; return true;
+  }
+
+// IMPLEMENTATIONS ARE FOUND IN
+#include "Poly.tcc"
+
+CORE_END_NAMESPACE
+#endif

Packages/Core/include/CORE/poly/Poly.tcc

@@ -0,0 +1,685 @@
+/******************************************************************
+ * Core Library Version 1.6, June 2003
+ * Copyright (c) 1995-2003 Exact Computation Project
+
+ File: Poly.tcc
+ Purpose: 
+	Template implementations of the functions
+	of the Polynomial<NT> class (found in Poly.h)
+
+ OVERVIEW:  
+	Each polynomial has a nominal "degree" (this
+		is an upper bound on the true degree, which
+		is determined by the first non-zero coefficient).
+	Coefficients are parametrized by some number type "NT".
+	Coefficients are stored in the "coeff" array of
+		length "degree + 1".  
+	IMPORTANT CONVENTION: the zero polynomial has degree -1
+		while nonzero constant polynomials have degree 0.
+	
+ Author: Chee Yap, Sylvain Pion and Vikram Sharma
+ Date:   May 28, 2002
+ Core Library
+ $Id$
+ ************************************** */ 
+
+template <class NT>
+  int Polynomial<NT>::COEFF_PER_LINE  = 3;           // pretty print parameters
+template <class NT>
+  const char* Polynomial<NT>::INDENT_SPACE ="   ";  // pretty print parameters
+
+// ==================================================
+// Polynomial Constructors
+// ==================================================
+
+template <class NT>
+  Polynomial<NT>::Polynomial(void) {
+    degree = -1; 	// this is the zero polynomial!
+    coeff = NULL;
+  }
+
+template <class NT>
+  Polynomial<NT>::Polynomial<NT>(int n) {
+    assert(n>= -1);
+    degree = n;
+    if (n == -1) return;	// return the zero polynomial!
+    if (n>=0) coeff = new NT[n+1];
+    coeff[0]=1;			// otherwise, return the unity polynomial
+    for (int i=1; i<=n; i++)
+      coeff[i]=0;
+  }
+
+template <class NT>
+  Polynomial<NT>::Polynomial<NT>(int n, NT * c) {
+    assert("array c has n+1 elements");
+    degree = n; 
+    if (n >= 0) {
+	coeff = new NT[n+1];
+	for (int i = 0; i <= n; i++) 
+	  coeff[i] = c[i];
+    }
+  }
+
+//
+// Constructor from a vector of NT's
+///////////////////////////////////////
+template <class NT>
+  Polynomial<NT>::Polynomial<NT>(const VecNT & vN) {
+    degree = vN.size()-1; 
+    if (degree >= 0) {
+	coeff = new NT[degree+1];
+	for (int i = 0; i <= degree; i++) 
+	  coeff[i] = vN[i];
+    }
+  }
+
+template <class NT>
+  Polynomial<NT>::Polynomial<NT>(const Polynomial<NT> & p) {
+    coeff = NULL;
+    *this = p;	// reduce to assignment operator=
+  }
+
+// Constructor from a Character string of coefficients
+///////////////////////////////////////
+template <class NT>
+  Polynomial<NT>::Polynomial<NT>(int n, const char * s[]) {
+    assert("array s has n+1 elements");
+    degree = n;
+    if (n >= 0) {
+        coeff = new NT[n+1];
+        for (int i = 0; i <= n; i++)
+          coeff[i] = s[i];
+    }
+  }
+
+template <class NT>
+  Polynomial<NT>::~Polynomial<NT>() {
+     if (degree >= 0)
+	delete[] coeff;
+  }
+
+  // assignment:
+template <class NT>
+  Polynomial<NT> & Polynomial<NT>::operator=(const Polynomial<NT>& p) {
+    if (this == &p) return *this;	// self-assignment
+    degree = p.getDegree();
+    delete[] coeff;
+    coeff = new NT[degree +1];
+    for (int i = 0; i <= degree; i++)
+	coeff[i] = p.coeff[i];
+    return *this;
+  }
+
+  // getTrueDegree
+template <class NT>
+  int Polynomial<NT>::getTrueDegree() const {
+    for (int i=degree; i>=0; i--)
+	if (coeff[i] != 0) return i;
+    return -1;	// Zero polynomial
+  }
+
+// ==================================================
+// polynomial arithmetic
+// ==================================================
+
+  // Expands the nominal degree to n;  
+  //	Returns n if nominal degree is changed to n
+  //	Else returns -2
+template <class NT>
+  int Polynomial<NT>::expand(int n) { 
+    if ((n <= degree)||(n < 0)) return -2;
+    int i;
+    NT * c = coeff;
+    coeff = new NT[n+1];
+    for (i = 0; i<= degree; i++) 
+      coeff[i] = c[i];
+    for (i = degree+1; i<=n; i++)
+      coeff[i] = 0;
+    delete[] c;
+    degree = n;
+    return n;
+  }
+
+  // contract() gets rid of leading zero coefficients
+  //	and returns the new (true) degree;
+  //	It returns -2 if this is a no-op
+template <class NT>
+  int Polynomial<NT>::contract() { 
+    int d = getTrueDegree();
+    if (d == degree) return (-2);  // nothing to do
+    else degree = d;
+    NT * c = coeff;
+    coeff = new NT[d+1];
+    for (int i = 0; i<= d; i++)
+	coeff[i] = c[i];
+    delete[] c;
+    return d;
+  }
+
+template <class NT>
+  Polynomial<NT> & Polynomial<NT>::operator+=(const Polynomial<NT>& p) { // +=
+    int d = p.getDegree();
+    if (d > degree) expand(d);
+    for (int i = 0; i<=d; i++)
+	  coeff[i] += p.coeff[i];
+    return *this;
+  }
+template <class NT>
+  Polynomial<NT> & Polynomial<NT>::operator-=(const Polynomial<NT>& p) { // -=
+    int d = p.getDegree();
+    if (d > degree) expand(d);
+    for (int i = 0; i<=d; i++)
+	  coeff[i] -= p.coeff[i];
+    return *this;
+  }
+
+  // SELF-MULTIPLICATION
+  // This is quadratic time multiplication!
+template <class NT>
+  Polynomial<NT> & Polynomial<NT>::operator*=(const Polynomial<NT>& p) { // *=
+    int d = degree + p.getDegree();
+    NT * c = new NT[d+1];
+    for (int i = 0; i<=d; i++)
+	c[i] = 0;
+
+    for (int i = 0; i<=p.getDegree(); i++)
+      for (int j = 0; j<=degree; j++) {
+	  c[i+j] += p.coeff[i] * coeff[j];
+      }
+    degree = d;
+    delete[] coeff;
+    coeff = c;
+    return *this;
+  }
+
+  // Multiply by a scalar
+template <class NT>
+  Polynomial<NT> & Polynomial<NT>::mulScalar( const NT & c) {
+    for (int i = 0; i<=degree ; i++)
+	coeff[i] *= c;
+    return *this;
+  }
+
+  // mulXpower: Multiply by X^i (COULD be a divide if i<0)
+template <class NT>
+  Polynomial<NT> & Polynomial<NT>::mulXpower(int s) {
+		// if s >= 0, then this is equivalent to 
+		// multiplying by X^s;  if s < 0, to dividing by X^s
+    if (s==0) return *this;
+    int d = s+getTrueDegree();
+    if (d < 0) {
+	degree = -1; delete[] coeff;
+        coeff = NULL;
+	return *this;
+    }
+    NT * c = new NT[d+1];
+    if (s>0) 
+	for (int j=0;  j <= d; j++) {
+	  if (j <= degree)
+		c[d-j] = coeff[d-s-j];
+	  else
+		c[d-j] = 0;
+	}
+    if (s<0) {
+	for (int j=0; j <= d; j++)
+	  c[d-j] = coeff[d-s-j];  // since s<0, (d-s-j) > (d-j)
+    }
+    delete[] coeff;
+    coeff = c; degree = d;
+    return *this;
+  }//mulXpower
+
+  // REDUCE STEP (helper for PSEUDO-REMAINDER function)
+  // Let THIS=(*this) be the current polynomial, and P be the input
+  //	argument for reduceStep.  Let R be returned polynomial.
+  //	R has the special form as a binomial,
+  //		R = C + M
+  //	where C is a constant and M a monomial (= coeff * power of X).
+  //	Moreover, THIS is transsformed to a new polynomial THAT which
+  //	is given by
+  // 		THAT =  (C * THIS) - M * P
+  //	MOREOVER: deg(THAT) < deg(THIS) unless deg(P)>deg(Q).  
+  //	Basically, C is a power of the leading coefficient of P.
+template <class NT>
+  Polynomial<NT> Polynomial<NT>::reduceStep (
+		Polynomial<NT>& p) {
+	p.contract(); contract();	// first contract both polynomials
+	Polynomial<NT> q(p);		// q is initially a copy of p
+					//	but is eventually M*P
+	int pDeg  = q.degree; int myDeg = degree;
+	if (pDeg == -1)  return *(new Polynomial());  // Zero Polynomial
+			// NOTE: pDeg=0 is really an error condition!
+	if (myDeg < pDeg) return *(new Polynomial(0));  // Unity Polynomial
+			// i.e., C=1, M=0.
+	// Now (myDeg >= pDeg).  Start to form the Return Polynomial R=C+M
+	Polynomial<NT> R(myDeg - pDeg + 1);  // deg(M)= myDeg - pDeg
+	q.mulXpower(myDeg - pDeg);  	 // q is going to become M*P
+
+	NT myLC = coeff[myDeg];	  // LC means "leading coefficient"
+	NT qLC = q.coeff[myDeg];  // p also has degree "myDeg" (qLC non-zero)
+	NT LC;
+	if ((myLC % qLC) == 0) { // myLC is divisible by qLC 
+	    LC = myLC / qLC;	 // exact division, LC not zero
+	    R.setCoeff(0, 1);  		  //  C = 1,
+	    if (LC != 1) 		  // Just a little optimization
+	        R.setCoeff(R.degree, LC); //  M = LC * X^(myDeg-pDeg)
+		q.mulScalar(LC); 	  //  q = M*P.
+	}
+	if ((qLC % myLC) == 0) {  // qLC is divisible by myLC 
+	    LC = qLC / myLC;	  // exact division, LC not zero
+	    if ((LC != 1) && (LC != -1)) { // IMPORTANT: unlike the previous
+				  // case, we need an additional condition
+				  // that LC != -1.  THis is because 
+				  // if (LC = -1), then we have qLC and
+				  // myLC are mutually divisible, and
+				  // we would be updating R twice! 
+		R.setCoeff(0, LC); 	   // C = LC, M = X^(myDeg-pDeg)
+		mulScalar(LC); 	   	   // And q = M*P.
+	    }
+	}
+	if (((qLC % myLC) != 0) && ((myLC % qLC)!= 0)) {  //
+	  // NT g = gcd(qLC, myLC); 	// This ASSUMES gcd is defined in NT !!
+	  NT g = 1;  			// In case no gcd is available
+	  if (g == 1) {
+	  	R.setCoeff(0, qLC);	  	// C = qLC
+	  	R.setCoeff(R.degree, myLC);	 // M = (myLC) * X^{myDeg-pDeg}
+	  	mulScalar(qLC);	 		// forming  C * THIS 
+	  	q.mulScalar(myLC);		// forming  M * P
+	  } else {
+	  	R.setCoeff(0, qLC/g);	  	// C = qLC/g
+	  	R.setCoeff(R.degree, myLC/g);	// M = (myLC/g) * X^{myDeg-pDeg}
+	  	mulScalar(qLC/g);	 	// forming  C * THIS 
+	  	q.mulScalar(myLC/g);		// forming  M * P
+	  }
+	} 
+	(*this) -= q;		// THAT is finally C*THIS - M*P
+	contract();		// Leading term is eliminated!
+	return R;		// Returns R = C + M
+  }// reduceStep
+
+  // PSEUDO-REMAINDER (and PSEUDO-QUOTIENT)
+  // Let THIS = (*this) be input polynomial and it is eventually
+  // transformed to THAT.  Let R be the returned polynomial,
+  // and P be the input polynomial.  Then we have
+  // 		THAT =  (C * THIS) - R * P
+  // where deg(THAT) < deg(P).
+  // Moreover, C is uniquely determined (we won't spell it out)
+  // except to note that
+  //	C divides D = (LC)^{deg(THIS)-deg(P)+1} 
+  //	where LC is the leading coefficient of P.
+  // In this program, C is stored as the static variable ccc_ (HACK)
+  // NOTE: Normally, Pseudo-Remainder is defined so that C = D
+
+template <class NT>
+  Polynomial<NT> Polynomial<NT>::pseudoRemainder (
+		Polynomial<NT>& p) {
+    contract(); p.contract();
+    Polynomial<NT> q(p);
+    if (q.degree == -1)  {
+	std::cout << "ERROR in Polynomial<NT>::pseudoRemainder :\n" <<
+		"    -- divide by zero polynomial" << std::endl;
+	return Polynomial(0);  // Unit Polynomial (arbitrary!)
+    }
+    if (q.degree > degree) {
+	return Polynomial(); // Zero Polynomial
+		// CHECK: THAT = 1 * THIS - 0 * P,  deg(THAT) < deg(P)
+    }
+    Polynomial<NT> R;  	// accumulating the return polynomial R
+    Polynomial<NT> tmpR;
+    int sgn = -1;
+//    ccc_ = *(new NT(1));  // ccc_ is a static NT variable (hack)
+    while (degree >= p.degree) {
+        q = p;
+	tmpR = reduceStep(q);		// tmpR = C + M
+        if (tmpR.coeff[0] < 0)
+          sgn = -sgn;
+	R.mulScalar(tmpR.coeff[0]);	// R *= C
+//	ccc_ *= tmpR.coeff[0];		// hack!!
+	tmpR.setCoeff(0,0);		// remove the zeroth coeff
+	R += tmpR.mulXpower(-1);	// R += M   (so, R -> C*R + M)
+    }
+
+    if (sgn<0) *this = -*this;
+    return R;	// R is the pseudo-quotient
+  }//pseudoRemainder
+
+template <class NT>
+  Polynomial<NT> & Polynomial<NT>::operator-() {	// unary minus
+    for (int i=0; i<=degree; i++)
+	coeff[i] *= -1;
+    return *this;
+  }
+template <class NT>
+  Polynomial<NT> & Polynomial<NT>::power(unsigned int n) {	// self-power
+    if (n == 0) {
+	degree = 0;
+	delete [] coeff;
+	coeff = new NT[1]; coeff[0] = 1;
+    } else {
+      Polynomial<NT> p = *this;
+      for (unsigned int i=1; i<n; i++)
+	*this *= p;
+    }
+    return *this;
+  }
+
+  // evaluation of Expr value
+template <class NT>
+  Expr Polynomial<NT>::eval(const Expr& e) const {		// evaluation
+    if (degree == -1) return Expr(0);
+    if (degree == 0) return Expr(coeff[0]);
+    Expr val(0);
+    for (int i=degree; i>=0; i--) {
+	val *= e;
+	val += coeff[i];
+    }
+    return val;
+  }//eval
+
+  // evaluation of BigFloat value
+template <class NT>
+  BigFloat Polynomial<NT>::eval(const BigFloat& f) const {		// evaluation
+    if (degree == -1) return BigFloat(0);
+    if (degree == 0) return BigFloat(coeff[0]);
+    BigFloat val(0);
+    for (int i=degree; i>=0; i--) {
+	val *= f;
+	val.makeExact();	// in future, this should NOT be necessary
+	val += coeff[i];
+	val.makeExact();	// in future, this should NOT be necessary
+    }
+    return val;
+  }//eval
+
+  //============================================================
+  // Bounds
+  //============================================================
+
+  // Cauchy Upper Bound on Roots
+template < class NT >
+  BigFloat Polynomial<NT>::CauchyUpperBound() const {
+    if (zeroP(*this))
+      return BigFloat(0);
+    Expr mx = 0;
+    int deg = getTrueDegree();
+    for (int i = 0; i < deg; ++i) {
+      mx = core_max(mx, Expr(core_abs(coeff[i])));
+    }
+    mx /= core_abs(coeff[deg]);
+    mx.approx(CORE_INFTY, 2); 
+    	// get an absolute approximate value with error < 1/4
+    return (mx.BigFloatValue().makeExact() + 2);
+}
+  // Cauchy Lower Bound on Roots
+template < class NT >
+  BigFloat Polynomial<NT>::CauchyLowerBound() const {
+    if (zeroP(*this))
+      return BigFloat(0);
+    Expr mx = 0;
+    int deg = getTrueDegree();
+    for (int i = 1; i <= deg; ++i) {
+      mx = core_max(mx, Expr(core_abs(coeff[i])));
+    }
+    mx = core_abs(coeff[0])/(core_abs(coeff[0]) + mx) ;
+    mx.approx(2, CORE_INFTY); 
+    	// get an relative approximate value with error < 1/4
+    return (mx.BigFloatValue().makeExact().div2());
+}
+
+  // Separation bound for polynomials that may have multiple roots
+  // We use the Rump-Schwartz bound
+template < class NT >
+  BigFloat Polynomial<NT>::sepBound() const {
+    BigInt d;
+    BigFloat e;
+    int deg = getTrueDegree();
+    
+    CORE::power(d, BigInt(deg), ((deg)+5)/2);
+    e = CORE::power(height()+1, deg);
+    return 1/e*2*d;
+  }
+
+  // height function
+template < class NT >
+  BigFloat Polynomial<NT>::height() const {
+    if (zeroP(*this))
+	return NT(0);
+    int deg = getTrueDegree();
+    NT ht = 0;
+    for (int i = 0; i< deg; i++) 
+	ht = std::max(ht, core_abs(coeff[i]) );
+    return BigFloat(ht); 
+  }
+
+  // length function
+template < class NT >
+  BigFloat Polynomial<NT>::length() const {
+   if (zeroP(*this))
+	return NT(0);
+    int deg = getTrueDegree();
+    NT length = 0;
+    for (int i = 0; i< deg; i++) 
+	length += core_abs(coeff[i]*coeff[i]);
+    return sqrt(BigFloat(length));
+  }
+
+  //============================================================
+  // differentiation
+  //============================================================
+
+template <class NT>
+  Polynomial<NT> & Polynomial<NT>::differentiate() {	// self-differentiation
+    if (degree >= 0) {
+      NT * c = new NT[degree];
+      for (int i=1; i<=degree; i++) 
+	c[i-1] = coeff[i] * i;
+      degree--;
+      delete[] coeff; 
+      coeff = c;
+    }
+    return *this;
+  }// differentiation
+
+  // multi-differentiate
+template <class NT>
+  Polynomial<NT> & Polynomial<NT>::differentiate(int n) {
+    assert(n >= 0);
+    for (int i=1; i<=n; i++)
+	this->differentiate();
+    return *this;
+  } // multi-differentiate
+
+
+// ==================================================
+// Square-free and Primitive Parts
+// ==================================================
+
+// squareFree not yet implemented
+
+
+// Primitive Part:  *this is changed
+//
+template <class NT>
+  Polynomial<NT> & Polynomial<NT>::primPart(){
+  	// NOTE: GCD must be provided by NT
+    int d = getTrueDegree();
+    assert (d >= 0);
+    if (d == 0) {
+	if (coeff[0] > 0)
+	   coeff[0] = 1;
+	else
+	   coeff[0] = -1;
+	return *this;
+    }
+
+    NT g = coeff[d];
+    NT c;
+    for (int i=d-1; i>=0; i--) {
+    	c = coeff[i];
+        g = gcd(g, c);
+    }
+    if (g==1) return *this;
+
+    for (int i=0; i<=d; i++){
+    	coeff[i] =  coeff[i]/g;
+    }
+    return *this;
+  }// primPart
+
+
+// ==================================================
+// Useful member functions
+// ==================================================
+
+  // reverse:
+  // 	reverses the list of coefficients
+template <class NT>
+  void Polynomial<NT>::reverse() { 
+    NT tmp;
+    for (int i=0; i<= degree/2; i++) {
+  	tmp = coeff[i];
+  	coeff[i] =   coeff[degree-i];
+  	coeff[degree-i] = tmp;
+    }
+  }//reverse
+
+  // dump: prints out polynomial
+template <class NT>
+  void Polynomial<NT>::dump() const {
+    std::cerr << "Polynomial: deg = " << degree << std::endl;
+    int termsInLine = 1;
+    int i=0;
+    for (; i<= getTrueDegree();  ++i) 
+	if (coeff[i] != 0) break;
+    
+    std::cerr << Polynomial<NT>::INDENT_SPACE << "(" << coeff[i] << ")";
+    if (i>0)
+      std::cerr << " * x ^ " << i;
+
+
+    for (i++ ; i<= getTrueDegree(); ++i) {
+      if (coeff[i] == 0) continue;
+      termsInLine++;
+      if (termsInLine % Polynomial<NT>::COEFF_PER_LINE == 0)
+	  std::cerr << std::endl << Polynomial<NT>::INDENT_SPACE ;
+      std::cerr << " + (" << coeff[i] << ")";
+      std::cerr << " * x ^ " << i;
+    }
+    std::cerr << std::endl;
+  }//dump
+
+  // dump: with optional message string
+template <class NT>
+  void Polynomial<NT>::dump(std::string m1) const{
+    std::cerr << m1;
+    dump();
+  }
+
+  // Dump of Maple Code for Polynomial
+template <class NT>
+  void Polynomial<NT>::mapleDump() const {
+    if (zeroP(*this)) {
+      std::cerr << 0 << std::endl;
+      return;
+    }
+    std::cerr << coeff[0];
+    for (int i = 1; i<= getTrueDegree(); ++i) {
+      std::cerr << " + (" << coeff[i] << ")";
+      std::cerr << "*x^" << i;
+    }
+    std::cerr << std::endl;
+  }//mapleDump
+
+// ==================================================
+// Useful friend functions for Polynomial<NT> class
+// ==================================================
+
+  // friend differentiation
+template <class NT>
+  Polynomial<NT> differentiate(const Polynomial<NT> & p) {	  // differentiate
+    Polynomial<NT> q(p);
+    return q.differentiate();
+  }
+
+  // friend multi-differentiation
+template <class NT>
+  Polynomial<NT> differentiate(const Polynomial<NT> & p, int n){//multi-differentiate
+    Polynomial<NT> q(p);
+    assert(n >= 0);
+    for (int i=1; i<=n; i++)
+	q.differentiate();
+    return q;
+  }
+
+  // friend equality comparison
+template <class NT>
+  bool operator==(const Polynomial<NT>& p, const Polynomial<NT>& q) {	// ==
+    int d, D;
+    Polynomial<NT> P(p); P.contract();
+    Polynomial<NT> Q(q); Q.contract();
+    if (P.degree < Q.degree) {
+	d = P.degree; D = Q.degree;
+	for (int i = d+1; i<=D; i++)
+	  if (Q.coeff[i] != 0) return false;	// return false
+    } else {
+	D = P.degree; d = Q.degree; 
+	for (int i = d+1; i<=D; i++)
+	  if (P.coeff[i] != 0) return false;	// return false
+    }
+    for (int i = 0; i <= d; i++)
+	if (P.coeff[i] != Q.coeff[i]) return false;	// return false
+    return true; 	// return true
+  }
+
+  // friend non-equality comparison
+template <class NT>
+  bool operator!=(const Polynomial<NT>& p, const Polynomial<NT>& q) {	// !=
+    return (!(p == q));
+  }
+
+  // stream i/o
+template <class NT>
+  std::ostream& operator<<(std::ostream& o, const Polynomial<NT>& p) {
+    o <<   "Polynomial<NT> ( deg = " << p.degree ;
+    if (p.degree >= 0) {
+      o << "," << std::endl;
+      o << ">  coeff c0,c1,... = " << p.coeff[0];
+      for (int i=1; i<= p.degree; i++)
+	o << ", " <<  p.coeff[i] ;
+    }
+    o << ")" << std::endl;
+    return o;
+  }
+
+  // fragile version...
+template <class NT>
+  std::istream& operator>>(std::istream& is, Polynomial<NT>& p) {
+    is >> p.degree; 
+    p.coeff = new NT[p.degree+1];
+    for (int i=0; i<= p.degree; i++)
+      is >> p.coeff[i];
+    return is;
+  }
+
+// ==================================================
+// Simple test of poly
+// ==================================================
+
+template <class NT>
+bool testPoly() {
+  int c[] = {1, 2, 3};
+  Polynomial<NT> p(2, c);
+  std::cout << p;
+
+  Polynomial<NT> zeroP;
+  std::cout << "zeroP  : " << zeroP << std::endl;
+
+  Polynomial<NT> P5(5);
+  std::cout << "Poly 5 : " << P5 << std::endl;
+
+  return 0;
+}
+
+// ==================================================
+// End of Polynomial<NT>
+// ==================================================

Packages/Core/include/CORE/poly/Sturm.h

@@ -0,0 +1,685 @@
+/******************************************************************
+ * Core Library Version 1.6, June 2003
+ * Copyright (c) 1995-2003 Exact Computation Project
+
+   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.
+> 
+>    TODO LIST and Potential Bugs:
+> 	(1) Split an isolating interval to give definite sign
+> 	(2) Should have a test for square-free polynomials
+> 	(3) The copy constructor seems to fail sometimes?
+
+   Author:  Chee Yap and Sylvain Pion, Vikram Sharma
+   Date:    July 20, 2002
+
+   Core Library  Version 1.5
+   $Id$
+ ************************************** */
+
+#include "../BigFloat.h"
+#include "Poly.h"
+
+CORE_BEGIN_NAMESPACE
+
+// ==================================================
+// Sturm Class
+// ==================================================
+
+template < class NT >
+class Sturm {
+public:
+  int len;			// length of array
+  Polynomial<NT> * seq;		// array of polynomials of length "len"
+  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
+  Sturm() : len(0), NEWTON_DIV_BY_ZERO(false){		// null constructor
+  }
+
+  Sturm(Polynomial<NT> pp) : NEWTON_DIV_BY_ZERO(false) { // constructor from polynomial
+    int ell = pp.getTrueDegree();
+    seq = new Polynomial<NT> [ell+1];
+    seq[0] = pp;	// 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
+    seq[1] = differentiate(pp);
+    int i;
+    for (i=2; i <= ell; i++) {
+        seq[i] = seq[i-2];
+	seq[i].pseudoRemainder(seq[i-1]);
+	if (zeroP(seq[i])) break;
+	seq[i].primPart();
+    }
+    len = i;
+  }
+
+  // copy constructor
+  Sturm(const Sturm&s) : len(s.len), NEWTON_DIV_BY_ZERO(s.NEWTON_DIV_BY_ZERO) {
+    seq = new Polynomial<NT> [len];
+    for (int i=0; i<len; i++)
+      seq[i] = s.seq[i];
+  }
+
+  ~Sturm() {
+    delete[] seq;
+  }
+
+  // METHODS
+
+  // dump functions
+  void dump(std::string msg) const {
+    std::cerr << msg << std::endl;
+    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 the first polynomial at x
+  //	PRE-CONDITION: sx != 0
+  int signVariations(const BigFloat & x, int sx) const
+  {
+    assert(sx != 0);
+    int cnt = 0;
+    int last_sign = sx;
+    for (int i=0; i<len; i++)
+    {
+      int sgn = seq[i].eval(x).sign();
+      if (sgn*last_sign < 0)
+        cnt++;
+      if (sgn != 0)
+        last_sign = sgn;
+    }
+    return cnt;
+  }
+
+  // signVariations(x)
+  //	--the first polynomial eval is not yet done
+  int signVariations(const BigFloat & x) const
+  {
+	int signx = seq[0].eval(x).sign();
+	if (signx == 0)
+		return (-1);	// error!
+	return signVariations(x, signx);
+  }//signVariations(x)
+
+  // signVariation at +Infinity
+  int signVariationsAtPosInfty() const
+  {
+    int cnt = 0;
+    int last_sign = seq[0].coeff[seq[0].getTrueDegree()].sign();
+    assert(last_sign != 0);
+    for (int i=1; i<len; i++)
+    {
+      int sgn = seq[i].coeff[seq[i].getTrueDegree()].sign();
+      if (sgn*last_sign < 0)
+        cnt++;
+      if (sgn != 0)
+        last_sign = sgn;
+    }
+    return cnt;
+  }
+
+  // signVariation at -Infinity
+  int signVariationsAtNegInfty() const
+  {
+    int cnt = 0;
+    int last_sign = seq[0].coeff[seq[0].getTrueDegree()].sign();
+    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 * seq[i].coeff[seq[i].getTrueDegree()].sign();
+      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
+  int numberOfRoots(const BigFloat &x, const BigFloat &y) const
+  {
+    // assert(x <= y);   // we allow x=y (probabably never used)
+    int signx = seq[0].eval(x).sign();
+    int signy = seq[0].eval(y).sign();
+    // 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())/2;
+    BigFloat newx, newy;
+    if (signx == 0) newx = x - sep; else newx = x;
+    if (signy == 0) newy = y + sep;  else newy = y;
+    return (signVariations(newx, seq[0].eval(newx).sign())
+    		- signVariations(newy, seq[0].eval(newy).sign()) );
+  }//numberOfRoots
+
+  // numberOfRoots above x:
+  //   NOTE: should do a special case for x=0
+  ///////////////////////////////////////////
+  int numberOfRootsAbove(const BigFloat &x) const
+  {
+    int signx = seq[0].eval(x).sign();
+    if (signx != 0)
+       return signVariations(x, signx) - signVariationsAtPosInfty();
+    BigFloat newx = x - (seq[0].sepBound())/2;
+    return signVariations(newx, seq[0].eval(newx).sign())
+		- signVariationsAtPosInfty();
+  }
+
+  // numberOfRoots below x:
+  //   NOTE: should do a special case for x=0
+  ///////////////////////////////////////////
+  int numberOfRootsBelow(const BigFloat &x) const
+  {
+    int signx = seq[0].eval(x).sign();
+    if (signx != 0) 
+       return signVariationsAtNegInfty() - signVariations(x, signx);
+    BigFloat newx = x + (seq[0].sepBound())/2;
+    return signVariationsAtNegInfty() 
+		- signVariations(newx, seq[0].eval(newx).sign());
+  }
+
+
+  // isolateRoots(x, y, v) 
+  ///   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
+   */
+  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)) { 	// if 0 is inside our interval
+	   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
+      		v.push_back(std::make_pair(x, y));
+    } else {
+      BigFloat mid = (x+y).div2();
+      isolateRoots(x, mid, v);
+      isolateRoots(mid, y, v);
+    }
+  }
+
+  // isolateRoots(v)
+  ///   isolates all roots and returns them in v
+  /**  	v is a vector of isolated intervals
+   */
+  void isolateRoots(BFVecInterval &v) const
+  {
+    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
+  BFInterval isolateRoot(int i) const
+  {
+    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) {
+      i += n+1;
+      if (i <= 0) return BFInterval(1,0); // ERROR CONDITION
+    }
+    if (n < i) return BFInterval(1,0);  // ERROR CONDITION INDICATED
+    if (n == 1) return BFInterval(x, y);
+    BigFloat m = (x+y).div2();
+    n = numberOfRoots(x, m);
+    if (n < i)
+      return isolateRoot(i-n, m, y);
+    else
+      return isolateRoot(i, x, m);
+  }
+
+  // same as isolateRoot(i).
+  BFInterval diamond(int i) const
+  {
+    return isolateRoot(i);
+  }
+
+  // First root above
+  BFInterval firstRootAbove(const BigFloat &e) const
+  {
+    return isolateRoot(1, e, seq[0].CauchyUpperBound());
+  }
+
+  // Main root (i.e., first root above 0)
+  BFInterval mainRoot() const
+  {
+    return isolateRoot(1, 0, seq[0].CauchyUpperBound());
+  }
+  
+  // First root below
+  BFInterval firstRootBelow(const BigFloat &e) const
+  {
+    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].eval(J.first) * seq[0].eval(J.second) < 0
+
+     BigFloat width = I.second - I.first;
+
+     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].eval(midpoint);
+       if (m == 0) {
+	       J.first = J.second = midpoint;
+	       return J;
+	}
+       if (seq[0].eval(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 experimental version of "refineAllRoots"
+  // 	based on Newton iteration
+  //
+  void newtonRefineAllRoots( 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 = 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
+
+
+  /////////////////////////////////////////////////////////////////
+  //  DIAGNOSTIC TOOLS 
+  /////////////////////////////////////////////////////////////////
+  // Polynomial tester:	P is polynomial to be tested
+  // 			prec is the bit precision for root isolation
+  // 			n is the number of roots predicted
+  
+  static void testSturm(PolyNT &P, int prec, int n = -1)
+  {
+    Sturm<NT> SS (P);
+    BFVecInterval v;
+    SS.refineAllRoots(v, prec);
+    cout << "   Number of roots is " << v.size();
+    if ((n >= 0) & (v.size() == (unsigned)n))
+    		cout << " (CORRECT!)" << endl;
+    else
+        	cout << " (ERROR!) " << endl;
+    int i = 0;
+    for (BFVecInterval::const_iterator it = v.begin();
+         	it != v.end(); ++it) {
+      	cout << ++i << "th Root is in ["
+		<< it->first << " ; " << it->second << "]" << endl;
+    }
+  }// testSturm
+
+  // testNewtonSturm( Poly, aprec, n)
+  // 	will run the Newton-Sturm refinement to isolate the roots of Poly
+  // 		until absolute precision aprec.
+  //	n is the predicated number of roots
+  //		(will print an error message if n is wrong)
+  static void testNewtonSturm(PolyNT &P, int prec, int n = -1)
+  {
+    Sturm<NT> SS (P);
+    BFVecInterval v;
+    SS.newtonRefineAllRoots(v, prec);
+    cout << "   Number of roots is " << v.size();
+    if ((n >= 0) & (v.size() == (unsigned)n))
+    		cout << " (CORRECT!)" << endl;
+    else
+        	cout << " (ERROR!) " << endl;
+
+    int i = 0;
+    for (BFVecInterval::const_iterator it = v.begin();
+         	it != v.end(); ++it) {
+      	cout << ++i << "th Root is in ["
+		<< it->first << " ; " << it->second << "]" << endl;
+    }
+  }// testNewtonSturm
+
+  // Newton(prec, bf):
+  // 	iterate until epsilon step size (i.e., epsilon < 2^{-prec})
+  //	starting from initial bf
+  //
+  BigFloat newton(long prec, const BigFloat& bf) const {
+	  // usually, prec is positive
+    int count = 0;
+    BigFloat val = bf;
+    extLong uMSB;
+    BigFloat del;
+    // newton iteration
+    do { 
+      del = seq[0].eval(val)/seq[1].eval(val);
+      del.makeExact();
+      uMSB = del.uMSB();
+      val -= del;       
+      val.makeExact();
+      count ++;
+    } while ((uMSB >= -prec) && (count < N_STOP_ITER)) ;
+			// N_STOP_ITER to prevent infinite loops
+    if (count == N_STOP_ITER) {
+	    core_error("Newton Iteration reach limit",
+			    __FILE__, __LINE__, true);
+    }
+    return val;
+  }
+
+  // v = newtonIterN(n, bf, del)
+  //    v is the root after n iterations of Newton
+  // 		starting from initial value of bf
+  // 	Return by reference, "del" (difference between returned val and value
+  // 		in the previous Newton iteration)
+  BigFloat newtonIterN(long n, const BigFloat& bf, BigFloat& del) {
+    int count = 0;
+    BigFloat val = bf;
+
+    extLong uMSB = val.uMSB();
+    // newton iteration
+    for (int i=0; i<n; i++) { 
+      BigFloat ff = seq[1].eval(val).makeExact();
+      if (ff == 0) {
+	 NEWTON_DIV_BY_ZERO = true;
+	 del = 0;
+	 core_error("Zero divisor in Newton Iteration", __FILE__, __LINE__, false);
+	 return 0;
+      }
+      BigFloat f = seq[0].eval(val).makeExact();
+      if (f == 0) {
+		 NEWTON_DIV_BY_ZERO = false;
+		 del = 0;    // indicates that we have reached the exact root
+		 return val; // val is the exact root, before the last iteration
+      }
+      del = (f/ff).makeExact();
+      uMSB = del.uMSB();
+      val -= del;
+      val.makeExact();
+      count ++;
+    }
+    return val;
+  }
+
+
+  // v = newtonIterE(prec, bf, del)
+  //
+  // 	return the value v which is obtained by Newton iteration until del.uMSB < -prec
+  // 		starting from initial value of bf
+  // 	Return by reference "del" (difference between returned val and value
+  // 		in the previous Newton iteration)
+
+  BigFloat newtonIterE(int prec, const BigFloat& bf, BigFloat& del) const{
+	// usually, prec is positive
+    int count = 0;
+    BigFloat val = bf;
+
+    do { 
+      val = newtonIterN(1, val, del);
+      count++;
+    } while ((del != 0) && ((del.uMSB() >= -prec) && (count < N_STOP_ITER))) ;
+				// N_STOP_ITER to prevent infinite loops
+    return val;
+  }
+
+  
+  // Smale's bound
+  //   smalesBound(p, q) where q is the derivative of p
+  //   returns Smale's original bound for convergence in range space
+  
+  BigFloat smaleBound(const Polynomial<NT> & p, const Polynomial<NT> & q) const{
+    int deg = p.getTrueDegree();
+    BigFloat rho = 1/(1 + pow(q.length(), deg)
+		    	  *pow(p.length() + 1, deg-1) );
+    return rho/13;
+  }
+
+  // Yap's bound
+  //   yapsBound(p) returns the bound on size of isolating interval
+  //   which will guarantee being in Newton zone
+  
+  BigFloat yapsBound(const Polynomial<NT> & p) const{
+    int deg = p.getTrueDegree();
+    return  1/(1 + pow(BigFloat(deg), 3*deg+9)
+	    	       *pow(BigFloat(2+p.height()),6*deg));
+  }
+
+  // newtonRefine(I, a) assumes I is an isolating interval for a root x^*,
+  //     and will return an approximate root x which is guaranteed
+  //     to be in the Newton basin of the root, and such that
+  //	|x-x^*| < 2^{-a}.
+
+  BFInterval newtonRefine(const BFInterval I, int aprec) {
+    BFInterval J = I;
+    int leftSign = seq[0].eval(I.first).sign();
+    if (leftSign == 0) { J.first = J.second = I.first; return J;}
+
+    int rightSign = seq[0].eval(I.second).sign();
+    if (rightSign == 0) { J.first = J.second = I.second; return J;}
+
+    assert( leftSign * rightSign < 0 );
+
+    int steps  = 1; //The number of times Newton is called without checking 
+    	// whether the iterations are still in the interval
+    BigFloat x;
+    BigFloat smale = smaleBound(seq[0], seq[1]);
+    BigFloat yap = yapsBound(seq[0]);
+    int N = steps;
+
+    BigFloat oldval = 0.0, temp;
+    BigFloat del;
+
+    x = (J.second + J.first).div2();
+    // x.makeExact();
+    while ( //seq[0].eval(x).abs() > smale && (J.second - J.first) > yap && 
+	   del.uMSB() >= - aprec ){
+	oldval = x;
+
+	x = newtonIterN(N, x, del);
+	if (del == 0) {  // either reached exact root or Divide-by-zero.  It is
+			// the user's responsibility to check which is the case.
+		J.first = x; J.second = x; 
+		return J;	
+	}
+	del = core_abs(del);
+	if  ( 
+	     (J.first <= x && x <= J.second)
+		        &&
+	     ((J.first < x - del) || (J.second > x+ del)))
+	{		// we can get a better root:
+	  if (x-del > J.first) {
+		  int  lSign = seq[0].eval(x - del).sign();
+		  if (lSign == 0) {
+			  J.first = J.second = x-del; return J;
+		  }
+		  if (lSign == leftSign)
+			  J.first = x - del;
+		  else {
+			  J.second = x - del;
+			  x = (J.second + J.first).div2();
+		  }
+	  }
+	  if (x+del < J.second) {
+		  int  rSign = seq[0].eval(x + del).sign();
+		  if (rSign == 0) {
+			  J.first = J.second = x+del; return J;
+		  }
+		  if (rSign == rightSign)
+			  J.second = x+del;
+		  else {
+		 	J.first = x + del;
+			x = (J.second + J.first).div2();
+		  }
+	  }
+
+	  N ++;		// be more confident or aggressive
+
+	} else {
+          N --;
+
+	  if (N == 0) {   // do bisection in this case
+	    x = (J.second + J.first).div2();  // x.makeExact();
+       		if(seq[0].eval(J.first) * seq[0].eval(x) < 0){
+	 		J.second = x;  
+       		}else if(seq[0].eval(J.second) * seq[0].eval(x) < 0){
+	 		J.first = x;
+       		}else{
+	 		J.first = x;
+	 		J.second = x;
+			return J;
+       		}
+       		x = (J.second + J.first).div2();  // x.makeExact();
+				//Reset x to be the midpoint of the interval
+       		N = steps;  	// restart
+
+	  } else {
+	    // Reset x to the oldval before doing Newton in this iteration
+	    	x = oldval;
+          }
+        }
+    }//end of outer while
+
+    /* For future work. We should exploit the fact that we have reached
+       the bounds. Could possible speed up things for high precision 
+    if(del.uMSB() >= -aprec){// If any of the bounds are satisfied 
+                             //before reaching the required precision
+      x = newtonIterE(aprec, x, del);
+      J.first = x; 
+      J.second = J.first;
+    }
+    */
+
+    return(J);
+
+  }//End of newton refine
+
+};// Sturm class
+
+// ==================================================
+// Static initialization
+// ==================================================
+template <class NT>
+int Sturm<NT>:: N_STOP_ITER = 10000;	// stop IterE after this many loops
+  					// Reset this as needed
+CORE_END_NAMESPACE
+