curves.h

CGAL is a collaborative effort of several sites in Europe and Israel. The goal is to make the most i

/****************************************************************************
 * Core Library Version 1.7, August 2004
 * Copyright (c) 1995-2004 Exact Computation Project
 * All rights reserved.
 *
 * This file is part of CORE (http://cs.nyu.edu/exact/core/); you may
 * redistribute it under the terms of the Q Public License version 1.0.
 * See the file LICENSE.QPL distributed with CORE.
 *
 * Licensees holding a valid commercial license may use this file in
 * accordance with the commercial license agreement provided with the
 * software.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * File: Curves.h
 *
 * Description: 
 * 	Two templated classes are defined here:
 *		Curve and BiPoly
 *	These classes are parametrized by the number type
 *		(called NT) which represents the
 *		domain of the coefficients of the underlying
 *		polynomials.  Standard default is NT=BigInt, but
 *		we will allow NT=int, NT=BigRat, NT=BigFloat, NT=Expr.
 *	BiPoly represents the class of bivariate polynomials,
 *		i.e.,  each BiPoly object is an element of NT[X,Y].
 *		We store each BiPoly as a list of polynomials in X.
 *	Curve represents the class of plane curves whose equation
 *		is A(X,Y)=0, for some BiPoly A(X,Y).
 *	Features:
 *		--Constructor from strings such as
 *			"3 x^2 + 7 xy^2 - 4 x + 13".
 *		--Basic plot functions
 *
 *	To Do:
 *	  --Dump should produce human readable strings like
 *	  	"3 x^2 + 7 xy^2 - 4 x + 13".
 *	  --String constructor generalizations:
 *	  	(1) allow one "=" sign (e.g., "3 x^2 = y^2 - xy")(DONE)
 *		(2) allow general parenthesis
 *		(3) allow X and Y (DONE)
 *	  --We should be able to read/write
 *	  	curve definitions from/to files
 *	  --Plot should be more efficient (use previous roots
 *	  	to help find the next roots, there should be
 *	  	a "plot structure" that is persistent)
 *	  --Plot should refine in both x- and y-increments.
 *	  --Plot should have some option to show the
 *	  	x- and y-axes, and to label some points.
 *	  --verticalIntersect(...) should be implemented using
 *	        Polynomial<BigFloat>, not Polynomial<Expr> for efficiency
 *	  --the plot parameters (eps,xmin,xmax,ymin,ymax) must be
 *	        made part of the Curve class (static members).
 *	        Incorporate the "setParams" method into class.
 *
 *  Author:  Vikram Sharma and Chee Yap
 *  Date:    April 12, 2004
 *
 * WWW URL: http://cs.nyu.edu/exact/
 * Email: [e-mail]
 *
 * $Source: /CVSROOT/CGAL/Packages/Core/include/CORE/poly/Curves.h,v $
 * $Revision: 1.2 $ $Date: 2004/11/14 12:00:18 $
 ***************************************************************************/


#ifndef CORE_CURVES_H
#define CORE_CURVES_H

#include <fstream>
#include <list>
#include "CORE/poly/Poly.h"

CORE_BEGIN_NAMESPACE

// ==================================================
// Curve Class
// ==================================================

//typedef BigInt NT;
//typedef Expr   NT;
//typedef Polynomial<NT>        PolyNT;
//typedef std::vector<Expr>	VecExpr;
//typedef std::vector<BigInt>	VecBigInt;
//typedef std::vector<NT>       VecNT;
//typedef std::vector<Polynomial<NT> >	VecPoly;

template <class NT>
class Monomial{
  //Helper class to store the coefficients for given x-deg and y-deg 
  //Used by string input routines
 public:
  NT coeff;
  int xdeg;
  int ydeg;

  Monomial(){
  }

  Monomial(NT& cf, int& dx, int& dy){
    coeff = cf;
    xdeg = dx;
    ydeg = dy;
  }

  void dump(){
    std::cout << coeff << "X^" << xdeg << " Y^" << ydeg;
  }
};


//Class of Bivariate polynomials
//	Viewed as a polynomial in Y with
//	coefficients which are polynomials in X
template <class NT>
class BiPoly{
 private:
  //The following are used in the constructor from strings.
  //For more details see the related constructor.
  void constructFromString(string& s, char myX='x', char myY='y');
  void constructX(int n, BiPoly<NT>& P);
  void constructY(int n, BiPoly<NT>& P);
  int getnumber(const char* c, int i, unsigned int len, BiPoly<NT> & P);
  bool isint(char c);
  int getint(const char* c, int i, unsigned int len, int & n);
  int matchparen(const char* cstr, int start);
  int getbasicterm(string s, BiPoly<NT> & P);
  int getterm(string s, BiPoly<NT> & P);

 public:
  int ydeg; //Y-degree of the polynomial
  std::vector<Polynomial<NT> > coeffX; //vector of (1+ydeg) polynomials in X
	// If ydeg = d, then the polynomial is F(X,Y) =
        //   (Y^d * coeffX[d]) + (Y^{d-1} * coeffX[d-1]) +...+ (coeffX[0]).

  ////////////////////////////////////////////////////////
  //Constructors
  ////////////////////////////////////////////////////////

  //BiPoly()
  BiPoly(); // zero polynomial

  //BiPoly(n)
  BiPoly(int n);// creates a BiPoly with nominal y-degree equal to n.

  //BiPoly(vp)
  BiPoly(std::vector<Polynomial<NT> > vp); // From vector of Polynomials

  //BiPoly(p, flag):
  //	if true, it converts polynomial p(X) into P(Y)
  // 	if false, it creates the polynomial Y-p(X)
  BiPoly(Polynomial<NT> p, bool flag=false);
  
  //BiPoly(deg, d[], C[]):
  //	Takes in a list of list of coefficients.
  //	Each cofficient list represents a polynomial in X
  //
  //  deg - ydeg of the bipoly
  //  d[] - array containing the degrees of each coefficient (i.e., X poly)
  //  C[] - list of coefficients, we use array d to select the
  //      coefficients
  BiPoly(int deg, int *d, NT *C);

  //BiPoly(String s, char myX, char myY)
  //  myX and myY are names of the two variables.
  //  Default values of myX and myY are 'x' and 'y'.
  //  The string s has the form "3 x^2 + 7 xy^2 - 4 x + 13"
  //
  //  For now, we assume no parentheses, * or =.
  
  BiPoly(const string& s, char myX='x', char myY='y');
  BiPoly(const char* s, char myX='x', char myY='y');

  // copy constructor
  BiPoly(const BiPoly<NT>&);

  //Destructor
  ~BiPoly();
  //Destructor helper
  void deleteCoeffX();


  ////////////////////////////////////////////////////////
  // METHODS
  ////////////////////////////////////////////////////////
  
  // filedump (msg, ofs, com, com2)
  // 	where msg, com, com2 are strings.
  // 	msg is an message and com, com2 are the strings
  // 	preceding each output line
  // 	(e.g., msg="BiVariate Polynomial"  and com=com2="# ")
  // This is called by the other dump functions
  void dump(std::ostream & os, std::string msg = "",
      std::string com="# ", std::string com2 = "# ") const;
  // dump(ofs, msg, com) -- dump to file
  //void dump(std::ofstream & ofs, std::string msg,
  //    std::string com="# ", std::string com2="# ") const;

  // dump(msg, com) -- dump to std output
  void dump(std::string msg="", std::string com="",
      std::string com2="") const;

  /*Cannot work with these two functions right now.
    BiPoly as per now can only handle BigInt and int since
    Expr cannot be handled by Polynomial class.*/
  
  // yPolynomial(x) 
  //   returns the polynomial (in Y) when we substitute X=x
  
  /* BiPoly<NT> yPolynomial(const Expr & x) {

    VecExpr vE;

    for (int i=0; i<= ydeg; i++) {
      vE.push_back(coeffX[i].eval(x));
    }
    
    return BiPoly<NT>(vE);
  }//yPolynomial
  */

  Polynomial<NT> yPolynomial(const NT & x);

  // Expr version of yPoly (temporary hack)
  Polynomial<Expr> yExprPolynomial(const Expr & x);

  // BF version of yPoly (temporary hack)
  Polynomial<BigFloat> yBFPolynomial(const BigFloat & x);

  // xPolynomial(y) 
  //   returns the polynomial (in X) when we substitute Y=y
  //   
  //   N.B. May need the
  //   		Polynomial<Expr> xExprPolynomial(Expr y)
  //   version too...
  //
  Polynomial<NT> xPolynomial(const NT & y) ;
  
  // getYdegree()
  int getYdegree() const;
  
  // getXdegree()
  int getXdegree();

  // getTrueYdegree
  int getTrueYdegree();

  //eval(x,y)
  Expr eval(Expr x, Expr y);//Evaluate the polynomial at (x,y)

  ////////////////////////////////////////////////////////
  // Polynomial arithmetic (these are all self-modifying)
  ////////////////////////////////////////////////////////
  
  // Expands the nominal y-degree to n;
  //	Returns n if nominal y-degree is changed to n
  //	Else returns -2

  int expand(int n);

  // contract() gets rid of leading zero polynomials