diff --git a/Packages/Core/include/CORE/geom2d/circle2d.h b/Packages/Core/include/CORE/geom2d/circle2d.h new file mode 100644 index 00000000000..93dd92b5499 --- /dev/null +++ b/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 + * Chen Li + * Zilin Du + * + * 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 diff --git a/Packages/Core/include/CORE/geom2d/line2d.h b/Packages/Core/include/CORE/geom2d/line2d.h new file mode 100644 index 00000000000..349ec1a6a90 --- /dev/null +++ b/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 + * Chen Li + * Zilin Du + * + * 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 diff --git a/Packages/Core/include/CORE/geom2d/point2d.h b/Packages/Core/include/CORE/geom2d/point2d.h new file mode 100644 index 00000000000..650e89064be --- /dev/null +++ b/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 + * Chen Li + * Zilin Du + * + * 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 diff --git a/Packages/Core/include/CORE/geom2d/segment2d.h b/Packages/Core/include/CORE/geom2d/segment2d.h new file mode 100644 index 00000000000..f583522c262 --- /dev/null +++ b/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 + * Chen Li + * Zilin Du + * + * 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 diff --git a/Packages/Core/include/CORE/poly/Poly.h b/Packages/Core/include/CORE/poly/Poly.h new file mode 100644 index 00000000000..c4bf114faa8 --- /dev/null +++ b/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 + +CORE_BEGIN_NAMESPACE + +class Expr; +// ================================================== +// Typedefs +// ================================================== + +//typedef std::vector VecExpr; +//typedef std::pair Interval; +//typedef std::vector VecInterval; +typedef std::pair BFInterval; + // NOTE: an error condition is indicated by + // the special interval (1, 0) +typedef std::vector BFVecInterval; + + +// ================================================== +// Polynomial Class +// ================================================== + +template +class Polynomial { +public: + typedef std::vector 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 & polyZero(); + static const Polynomial & 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::ccc_; + +// ================================================== +// Static Constants +// Does this belong here? +// ================================================== + +template < class NT > + CORE_INLINE + const Polynomial & Polynomial::polyZero(){ + static Polynomial zeroP; + return zeroP; + } + +template < class NT > + CORE_INLINE + const Polynomial & Polynomial::polyUnity() { + static NT c[] = {1}; + static Polynomial unityP(0, c); + return unityP; + } + +// ================================================== +// Useful functions for Polynomial class +// ================================================== + + // polynomial arithmetic: +template < class NT > + Polynomial operator+(const Polynomial&, const Polynomial&);// + +template < class NT > + Polynomial operator-(const Polynomial&, const Polynomial&);// - +template < class NT > + Polynomial operator*(const Polynomial&, const Polynomial&);// * +template < class NT > + Polynomial power(const Polynomial&, int n); // power +template < class NT > + Polynomial differentiate(const Polynomial&); // differentiate +template < class NT > + Polynomial differentiate(const Polynomial&, int n); // multi-differ. + + // comparisons +template < class NT > + bool operator==(const Polynomial&, const Polynomial&); // == +template < class NT > + bool operator!=(const Polynomial&, const Polynomial&); // != +template < class NT > + bool zeroP(const Polynomial &); // =Zero Poly? +template < class NT > + bool unitP(const Polynomial &); // =Unit Poly? + + // stream i/o +template < class NT > + std::ostream& operator<<(std::ostream&, const Polynomial&); +template < class NT > + std::istream& operator>>(std::istream&, Polynomial&); + +// ================================================== +// Inline Functions +// ================================================== + + // friend polynomial arithmetic: +template < class NT > + CORE_INLINE + Polynomial operator+(const Polynomial& p1, + const Polynomial& p2){ // + + return Polynomial(p1) += p2; + } +template < class NT > + CORE_INLINE + Polynomial operator-(const Polynomial& p1, + const Polynomial& p2){ // - + return Polynomial(p1) -= p2; + } +template < class NT > + CORE_INLINE + Polynomial operator*(const Polynomial& p1, + const Polynomial& p2){ // * + return Polynomial (p1) *= p2; + } +template < class NT > + CORE_INLINE + Polynomial power(const Polynomial& p, int n){ // power + return Polynomial(p).power(n); + } + +// equal to zero poly? +template < class NT > + CORE_INLINE + bool zeroP(const Polynomial & p){ // =Zero Poly? + return (p.getTrueDegree()== -1); + } +template < class NT > + CORE_INLINE + bool unitP(const Polynomial & p){ // =Unit Poly? + int d = p.getTrueDegree(); + return ((d == 0) && p.coeff[0]==1 ); + } + +// get functions +template < class NT > + CORE_INLINE + int Polynomial::getDegree() const { + return degree; + } +// get TRUE leading coefficient +template < class NT > + CORE_INLINE + const NT & Polynomial::getLeadCoeff() const { + return getCoeff(getTrueDegree()); + } + +// get last non-zero coefficient +template < class NT > + CORE_INLINE + const NT & Polynomial::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::getCoeffs() { + return &coeff; + } +template < class NT > + CORE_INLINE + const NT& Polynomial::getCoeff(int i) const { + //if (i > degree) return NULL; + assert(i <= degree); + return coeff[i]; + } +// set functions +template < class NT > + CORE_INLINE + bool Polynomial::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 diff --git a/Packages/Core/include/CORE/poly/Poly.tcc b/Packages/Core/include/CORE/poly/Poly.tcc new file mode 100644 index 00000000000..b3ba43d2b03 --- /dev/null +++ b/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 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 + int Polynomial::COEFF_PER_LINE = 3; // pretty print parameters +template + const char* Polynomial::INDENT_SPACE =" "; // pretty print parameters + +// ================================================== +// Polynomial Constructors +// ================================================== + +template + Polynomial::Polynomial(void) { + degree = -1; // this is the zero polynomial! + coeff = NULL; + } + +template + Polynomial::Polynomial(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 + Polynomial::Polynomial(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 + Polynomial::Polynomial(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 + Polynomial::Polynomial(const Polynomial & p) { + coeff = NULL; + *this = p; // reduce to assignment operator= + } + +// Constructor from a Character string of coefficients +/////////////////////////////////////// +template + Polynomial::Polynomial(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 + Polynomial::~Polynomial() { + if (degree >= 0) + delete[] coeff; + } + + // assignment: +template + Polynomial & Polynomial::operator=(const Polynomial& 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 + int Polynomial::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 + int Polynomial::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 + int Polynomial::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 + Polynomial & Polynomial::operator+=(const Polynomial& 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 + Polynomial & Polynomial::operator-=(const Polynomial& 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 + Polynomial & Polynomial::operator*=(const Polynomial& 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 + Polynomial & Polynomial::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 + Polynomial & Polynomial::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 + Polynomial Polynomial::reduceStep ( + Polynomial& p) { + p.contract(); contract(); // first contract both polynomials + Polynomial 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 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 + Polynomial Polynomial::pseudoRemainder ( + Polynomial& p) { + contract(); p.contract(); + Polynomial q(p); + if (q.degree == -1) { + std::cout << "ERROR in Polynomial::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 R; // accumulating the return polynomial R + Polynomial 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 + Polynomial & Polynomial::operator-() { // unary minus + for (int i=0; i<=degree; i++) + coeff[i] *= -1; + return *this; + } +template + Polynomial & Polynomial::power(unsigned int n) { // self-power + if (n == 0) { + degree = 0; + delete [] coeff; + coeff = new NT[1]; coeff[0] = 1; + } else { + Polynomial p = *this; + for (unsigned int i=1; i + Expr Polynomial::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 + BigFloat Polynomial::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::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::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::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::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::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 + Polynomial & Polynomial::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 + Polynomial & Polynomial::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 + Polynomial & Polynomial::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 + void Polynomial::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 + void Polynomial::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::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::COEFF_PER_LINE == 0) + std::cerr << std::endl << Polynomial::INDENT_SPACE ; + std::cerr << " + (" << coeff[i] << ")"; + std::cerr << " * x ^ " << i; + } + std::cerr << std::endl; + }//dump + + // dump: with optional message string +template + void Polynomial::dump(std::string m1) const{ + std::cerr << m1; + dump(); + } + + // Dump of Maple Code for Polynomial +template + void Polynomial::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 class +// ================================================== + + // friend differentiation +template + Polynomial differentiate(const Polynomial & p) { // differentiate + Polynomial q(p); + return q.differentiate(); + } + + // friend multi-differentiation +template + Polynomial differentiate(const Polynomial & p, int n){//multi-differentiate + Polynomial q(p); + assert(n >= 0); + for (int i=1; i<=n; i++) + q.differentiate(); + return q; + } + + // friend equality comparison +template + bool operator==(const Polynomial& p, const Polynomial& q) { // == + int d, D; + Polynomial P(p); P.contract(); + Polynomial 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 + bool operator!=(const Polynomial& p, const Polynomial& q) { // != + return (!(p == q)); + } + + // stream i/o +template + std::ostream& operator<<(std::ostream& o, const Polynomial& p) { + o << "Polynomial ( 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 + std::istream& operator>>(std::istream& is, Polynomial& 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 +bool testPoly() { + int c[] = {1, 2, 3}; + Polynomial p(2, c); + std::cout << p; + + Polynomial zeroP; + std::cout << "zeroP : " << zeroP << std::endl; + + Polynomial P5(5); + std::cout << "Poly 5 : " << P5 << std::endl; + + return 0; +} + +// ================================================== +// End of Polynomial +// ================================================== diff --git a/Packages/Core/include/CORE/poly/Sturm.h b/Packages/Core/include/CORE/poly/Sturm.h new file mode 100644 index 00000000000..dd9bc5405d3 --- /dev/null +++ b/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 * 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 PolyNT; + // CONSTRUCTORS + Sturm() : len(0), NEWTON_DIV_BY_ZERO(false){ // null constructor + } + + Sturm(Polynomial pp) : NEWTON_DIV_BY_ZERO(false) { // constructor from polynomial + int ell = pp.getTrueDegree(); + seq = new Polynomial [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 [len]; + for (int i=0; i= 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 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 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= -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 & p, const Polynomial & 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 & 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 +int Sturm:: N_STOP_ITER = 10000; // stop IterE after this many loops + // Reset this as needed +CORE_END_NAMESPACE +