curves.tcc

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.tcc
 *
 *  Description: 
 *  	This file contains the implementations of
 *		functions defined by Curves.h
 *	It is included through Curves.h
 *  	Please see Curves.h for more details.
 *
 *  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.tcc,v $
 * $Revision: 1.2 $ $Date: 2004/11/14 12:00:18 $
 ***************************************************************************/


//CONSTRUCTORS FOR THE BIPOLY CLASS
  ////////////////////////////////////////////////////////
  //Constructors
  ////////////////////////////////////////////////////////

template <class NT>
BiPoly<NT>::BiPoly(){ // zero polynomial
    ydeg = -1;
  }
  
  //BiPoly(n)

template <class NT>
BiPoly<NT>::BiPoly(int n){// creates a BiPoly with nominal y-degree equal to n.

     // To support this constructor, you need functions
     // that are equivalent to "setCoeff(i, PolyNT q)" in the Poly Class
     // You also want "getCoeff(i)" functions.
     // E.g. BiPoly<NT> p(10);
     //      PolyNT q(3, {1,2,3,4});
     //      p.setCoeff(10, q);
     //      p.setCoeff(0, PolyNT());
     //
    ydeg = n;
    if (n<0) return; // coeffX not needed
    for(int i=0; i<= ydeg; i++){
      Polynomial<NT> temp;
      coeffX.push_back(temp);      
    }
  }

  //BiPoly(vp)
template <class NT>
BiPoly<NT>::BiPoly(std::vector<Polynomial<NT> > vp){ 
    // From vector of Polynomials
    ydeg = vp.size() - 1;
    if(ydeg >=0){
      coeffX = vp;
    }
  }
  
  //BiPoly(p, flag):
  //	if true, it converts polynomial p(X) into P(Y)
  // 	if false, it creates the polynomial Y-p(X)

template <class NT>
BiPoly<NT>::BiPoly(Polynomial<NT> p, bool flag){
    if (flag){
      ydeg = p.getTrueDegree();
      if(ydeg >=0){
        for(int i=0; i<=ydeg; i++){
	  Polynomial<NT> temp(0);
	  temp.setCoeff(0, p.getCoeffi(i));
	  coeffX.push_back(temp);	// does STL make a copy of temp?
        }//for
      }//if
    } else {
      ydeg = 1;
      coeffX.push_back(p);
      coeffX.push_back(Polynomial<NT>::polyUnity());
    }//else
  }//BiPoly(p)


  //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

template <class NT>
BiPoly<NT>::BiPoly(int deg, int *d, NT *C){

    ydeg = deg;
    Polynomial<NT> temp;
    int max=0;
    for(int i=0; i <=deg; i++)
      max = core_max(d[i],max);

    NT *c = new NT[max];

    for(int i=0; i<= deg; i++){
      for(int j=0; j <=d[i]; j++)
	c[j] = C[i+j];
      temp = Polynomial<NT>(d[i],c);
      coeffX.push_back(temp);      
    }
    delete[] c;
  }//BiPoly(deg,d[],C[])


//The BNF syntax is the following:-
//    [bipoly] -> [term]| [term] +/- [bipoly]
//    [term] -> [basic term]|[basic term] [term]|[basic term]*[term]
//    [basic term] -> [number]|'x'|'y'|[basic term]'^'[number]
//                    | '(' [bipoly]')'|'-'
//Unary minus is treated as a basic term
template <class NT>
BiPoly<NT>::BiPoly(const char * s, char myX, char myY){
	string ss(s);
	constructFromString(ss, myX, myY);
}
template <class NT>
BiPoly<NT>::BiPoly(const string & s, char myX, char myY){
	string ss(s);
	constructFromString(ss, myX, myY);
}
template <class NT>
void BiPoly<NT>::constructFromString(string & s, char myX, char myY){
  if((myX != 'x' || myX != 'X') && (myY != 'y' || myY != 'Y')){
    //Replace myX with 'x' and myY with 'y' in s.
    unsigned int loc = s.find(myX, 0);
    while(loc != string::npos){
      s.replace(loc,1,1,'x');
      loc = s.find(myX, loc+1);
    }
    loc = s.find(myY, 0);
    while(loc != string::npos){
      s.replace(loc,1,1,'y');
      loc = s.find(myY, loc+1);
    }
  }
  (*this) =  (*this).getbipoly(s);

}

//Constructor from another BiPoly
template <class NT>
BiPoly<NT>::BiPoly(const BiPoly<NT> &P) : ydeg(P.ydeg), coeffX(P.coeffX) {
}

  //Destructor
template <class NT>
void BiPoly<NT>::deleteCoeffX(){
  coeffX.clear();
}

template <class NT>
BiPoly<NT>::~BiPoly(){
  if (ydeg >= 0)
    deleteCoeffX();
}




  ////////////////////////////////////////////////////////
  // METHODS USED BY STRING CONSTRUCTOR ABOVE
  ////////////////////////////////////////////////////////


//Sets the input BiPoly to X^n
template <class NT>
void BiPoly<NT>::constructX(int n, BiPoly<NT>& P){
  assert(n>= -1);
  P.deleteCoeffX();//Clear the present coeffecients
  Polynomial<NT> q(n);//Nominal degree n
  q.setCoeff(n,NT(1)); 
  if (n>0) q.setCoeff(0,NT(0));
  P.coeffX.push_back(q);
  P.ydeg = 0;
}

//Sets the input BiPoly to Y^n
template <class NT>
void BiPoly<NT>::constructY(int n, BiPoly<NT>& P){
  P.ydeg = 0;
  P.deleteCoeffX();//Clear the present coeffecients
  NT cs[] = {1};
  Polynomial<NT> q(0, cs);
  P.coeffX.push_back(q);//P is the unit bipoly. Now we just multiply Y^n
  P.mulYpower(n);
}

//Returns in P the coeffecient starting from start
template <class NT>
int BiPoly<NT>::getnumber(const char* c, int start, unsigned int len,
			  BiPoly<NT> & P){
  int j=0;
  char *temp = new char[len];
  while(isint(c[j+start])){
    temp[j]=c[j+start];j++;
  }
  temp[j] = '\0';
  NT cf = NT(temp);
  P.deleteCoeffX();
  Polynomial<NT> q(0);
  q.setCoeff(0, cf);
  P.coeffX.push_back(q);
  P.ydeg = 0;
  delete[] temp;

  return (j-1+start);//Exactly the length of the number
}

template <class NT>
bool BiPoly<NT>::isint(char c){
  if(c == '0' || c == '1' || c == '2' || c == '3' || c == '4' ||
     c == '5' || c == '6' || c == '7' || c == '8' || c == '9')
    return true;
  else
    return false;
}


//Returns as integer the number starting from start in c
template <class NT>
int BiPoly<NT>::getint(const char* c, int start, unsigned int len,
			  int & n){
  int j=0;
  char *temp = new char[len];
  while(isint(c[j+start])){
    temp[j]=c[j+start];j++;
  }
  temp[j] = '\n';
  n = atoi(temp);
  delete[] temp;

  return (j-1+start);//Exactly the length of the number
}

//Given a string starting with an open parentheses returns the place
// which marks the end of the corresponding closing parentheses.
//Strings of the form (A).
template <class NT>
int BiPoly<NT>::matchparen(const char* cstr, int start){
  int count = 0;
  int j=start;
  
  do{
    if(cstr[j] == '('){
      count++;
    }
    if(cstr[j] == ')'){
      count--;
    }
    j++;      
  }while(count != 0 );//j is one more than the matching ')'
  
  return j-1;
}


template <class NT>
int BiPoly<NT>::getbasicterm(string s, BiPoly<NT> & P){
  const char * cstr = s.c_str();
  unsigned int len = s.length();
  int i=0;
  //BiPoly<NT> * temp = new BiPoly<NT>();

  if(isint(cstr[i])){
    i = getnumber(cstr, i, len, P);
  }else if(cstr[i] == 'x'||cstr[i] == 'X'){
    constructX(1, P);
  }else if(cstr[i] == 'y'||cstr[i] == 'Y'){
    constructY(1, P);
  }else if(cstr[i] =='('){
    int oldi = i;
    i = matchparen(cstr, i);
    string t = s.substr(oldi+1, i -oldi -1);
    P = getbipoly(t);
  }else if(cstr[i] == '-'){//Unary Minus
    P.ydeg = 0;
    P.deleteCoeffX();//Clear the present coeffecients
    NT cs[] = {-1};
    Polynomial<NT> q(0, cs);
    P.coeffX.push_back(q);
  }else{
    std::cout <<"ERROR IN PARSING BASIC TERM" << std::endl;
  }
  //i+1 points to the beginning of next syntaxtic object in the string.
 if(cstr[i+1] == '^'){
    int n;
    i = getint(cstr, i+2, len, n);
    P.pow(n);
  }
  return i;
}


template <class NT>
int BiPoly<NT>::getterm(string s, BiPoly<NT> & P){
  unsigned int len = s.length();
  if(len == 0){// Zero BiPoly
    P = BiPoly<NT>();
    return 0;
  }
  unsigned int ind, oind;
  const char* cstr =s.c_str();
  string t;
  //P will be used to accumulate the product of basic terms.
  ind = getbasicterm(s, P);
  while(ind != len-1 && cstr[ind + 1]!='+' && cstr[ind + 1]!='-' ){
    //Recursively get the basic terms till we reach the end or see
    // a '+' or '-' sign.
    if(cstr[ind + 1] == '*'){
      t = s.substr(ind + 2, len - ind -2);
      oind = ind + 2;
    }else{
      t = s.substr(ind + 1, len -ind -1);
      oind = ind + 1;
    }

    BiPoly<NT> R;
    ind = oind + getbasicterm(t, R);//Because the second term is the offset in
                                     //t
    P *= R;
  }

  return ind;
}

template <class NT>
BiPoly<NT> BiPoly<NT>::getbipoly(string s){

    //Remove white spaces from the string
    unsigned int cnt=s.find(' ',0);
    while(cnt != string::npos){
      s.erase(cnt, 1);
      cnt = s.find(' ', cnt);
    }

    unsigned int len = s.length();
    if(len <= 0){//Zero Polynomial
      return BiPoly<NT>();
    }

    //To handle the case when there is one '=' sign
    //Suppose s is of the form s1 = s2. Then we assign s to
    //s1 + (-1)(s2) and reset len
    unsigned int loc;
    if((loc=s.find('=',0)) != string::npos){
      s.replace(loc,1,1,'+');
      string s3 = "(-1)(";
      s.insert(loc+1, s3);
      len = s.length();
      s.insert(len, 1, ')');
    }
    len = s.length();

    const char *cstr = s.c_str();
    string t;
    BiPoly<NT> P;
    // P will be the polynomial in which we accumulate the
    //sum and difference of the different terms.
    unsigned int ind;

    if(cstr[0] == '-'){
      t = s.substr(1, len);
      ind = getterm(t,P) + 1;
      NT negone(-1);
      P.mulScalar(negone);
    }else{
      ind = getterm(s, P);
    }

    unsigned int oind =0;//the string between oind and ind is a term
    while(ind != len -1){
      BiPoly<NT> R;
      t = s.substr(ind + 2, len -ind -2);
      oind = ind;
      ind = oind + 2 + getterm(t, R);
      if(cstr[oind + 1] == '+')
	P += R;
      else if(cstr[oind + 1] == '-')
	P -= R;
      else
	std::cout << "ERROR IN PARSING BIPOLY! " << std::endl;
    }

    return P;
}

  ////////////////////////////////////////////////////////
  // METHODS
  ////////////////////////////////////////////////////////
  
  // filedump (msg, ofs, com, com2)
  // 	where msg and com are strings.
  // 	msg is an message and com is the character preceding each line
  // 	(e.g., msg=""  and com=com2="# ")
  // This is called by the other dump functions
template <class NT>
void BiPoly<NT>::dump(std::ostream & os, std::string msg,
    std::string com, std::string com2) const {
    if (msg != "")
      os << msg << std::endl;
    if(ydeg == -1) {
      os << com << " Zero Polynomial" << std::endl;
      return;
    }
    bool first = true;
    os << com;
    for (int i=0; i <= ydeg; i++){
      if (!zeroP(coeffX[i])){
	if (i % 3 == 0) os << std::endl << com2 ;  // output 3 coefficients per line
	if (first) first = false;
	else os << " + ";
	os << "[";
	coeffX[i].filedump(os,"","",com2); // Note: first comment string is ""
	os << "]";
        if (i > 0) os <<" * y^"<< i ;
      }
    }
  }//dump

  // dump(ofs, msg, com, com2) -- dump to file
  //
/*
template <class NT>
void BiPoly<NT>::dump(std::ofstream & ofs, std::string msg,
              std::string com, std::string com2) const {
    dump(ofs, msg, com, com2);
  }
*/
  // dump(msg, com) -- dump to std output
template <class NT>
void BiPoly<NT>::dump(std::string msg, std::string com,
    std::string com2) const {
    dump(std::cout, msg, com, com2);
  }

  /* ***********************************************************
    We want to substitute X or Y with a general Expression or BigFloat
    into a BiPoly.  E.g., the following produces a Y-polynomial
    after replacing X by x:

  BiPoly<NT> yPolynomial(const Expr & x) {